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ABSTRACT 


Wall-Fluid  Instabilities  in  Compliant  Channels  Conveying  Developing  Flows 

Peter  Gavin  LaRose 

A  collapsed  lung  airway  or  flexible  tube  is  modelled  as  a  two-dimensional  channel  of  infinite 
length.  We  consider  the  linear  stability  of  this  system  conveying  a  developing  flow  profile,  which 
is  approximated  with  a  uniform  profile  and  perturbation  of  the  Blasius  profile  for  flow  over  a  flat 
plate.  Exact  and  asymptotic  solutions  are  found  for  the  uniform  profile.  For  the  perturbation 
profile  an  analytical  solution  is  found  for  long  waves  and  a  numerical  shooting  solution  for 
arbitrary  wave  lengths.  Results  for  the  uniform  and  perturbation  profiles  are  found  to  be  generally 
in  qualitative  agreement,  though  the  perturbation  profile  is  less  stable.  For  the  perturbation  profile 
we  find  a  long  wave  instability  which  is  absent  for  uniform  flow  and  hence  has  not  been  seen  in 
previous  channel  studies.  This  is  stabilized  by  increasing  the  elastance  of  the  wall,  but  other  wall 
properties  do  not  affect  the  critical  flow  speed  except  in  correction  terms.  Increasing  the  channel 
width  decreases  the  critical  flow  speed  for  the  instability,  but  increases  the  critical  flow  rate.  We 
hypothesize  that  this  is  related  to  the  tube  collapse  that  is  seen  in  this  type  of  system  prior  to  the 
appearance  of  an  oscillatory  (flutter)  instability.  The  finite  wave  length  (flutter)  instability  is 
destabilized  by  decreasing  wall  damping,  increasing  wall  inertia,  decreasing  wall  elastance  or 
flexural  rigidity,  and  decreasing  channel  width,  and  may  appear  independent  of  or  simultaneously 
with  the  long  wave  instability.  Comparisons  with  experimental  investigations  of  air  flow  in 
flexible  tubes  shows  that  the  theoretically  predicted  flutter  frequencies  are  in  good  agreement  with 
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those  observed  experimentally,  but  that  there  are  difficulties  in  comparing  the  predicted  critical 
flow  speeds,  due  to  the  tube  geometry  introduced  by  the  tube  collapse  that  precedes  flutter.  We 
investigate  the  possible  effect  of  this  geometry  using  finite  elements  software  to  model  flow  in 
the  collapsed  tube  cross-section.  Finally,  we  compare  the  experimentally  observed  critical  flow 
speeds  for  the  onset  of  collapse  with  our  theoretical  predictions  for  the  long  wave  instability, 
finding  qualitative  support  for  our  hypothesis  that  collapse  is  the  physical  manifestation  of  the 
long  wave  instability. 
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CHAPTER  1:  INTRODUCTION 


§1.1  Motivation 

In  this  dissertation  we  examine  the  stability  of  a  flexible,  fluid-conveying  two-dimensional 
channel.  This  is  a  model  for  the  many  systems  in  engineering  and  physiology  in  which  flows  are 
conveyed  by  tubes  that  are  to  some  degree  flexible.  In  engineering  applications  the  existence  of 
flow-induced  vibrations  or  larger  scale  motion,  e.g.  in  oil  pipelines,  may  attest  to  the  flexibility 
of  a  tube,  and  in  other  applications,  e.g.  physiological  tubes  such  as  the  lung  airways,  the 
flexibility  of  such  a  structure  is  more  obvious.  In  many  of  these  applications  the  tubes  are  short 
and  flow  speeds  large,  so  that  the  flows  do  not  become  fully  developed.  This  is  the  case  for  the 
application  in  which  we  are  primarily  interested,  the  lung  airways,  which  have  air  flow  through 
short,  thick-walled  tubes.  It  is  thought  that  wheezing  lung  sounds  are  symptomatic  of  an 
oscillatory  wall-fluid  instability  (flutter),  and  this  motivates  the  present  study. 

§1.2  Characteristics  of  flexible  tubes 

There  are  a  number  of  interrelated  steady  and  unsteady  behaviors  that  occur  in  flexible 
fluid-conveying  tubes  which  are  absent  in  rigid  systems.  If  the  pressure  difference  between  the 
exterior  and  interior  of  such  a  tube  is  sufficiently  large,  the  tube  may  collapse,  resulting  in  an 
obstruction  of  the  tube  cross-sectional  area.  With  increasing  pressure,  the  initially  circular  cross- 
section  becomes  oval,  then  dumbbell  shaped,  and  finally  reduces  to  a  pair  of  small  conduits 
connected  by  a  completely  closed  center  section  (Shapiro  1977).  For  very  flexible  tubes  the 
pressure  difference  required  to  precipitate  collapse  may  be  small  enough  that  easily  attainable 
changes  in  the  tube  environment,  or  the  pressure  drop  resulting  from  a  flow  in  the  tube  interior, 
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may  be  sufficient  to  cause  some  degree  of  reduction  in  the  tube  cross-sectional  area.  With  the 
occurrence  of  collapse,  the  flow  rate  through  the  tube  may  become  independent  of  the  conditions 
at  the  downstream  end  of  the  tube,  either  from  wave-speed  limitation  (Dawson  and  Elliot  1977; 
Shapiro  1977)  or  viscous  limitation  (Shapiro  1977).  While  the  flow  rate  remains  constant, 
however,  the  flow  speed  in  the  tube  may  still  increase  with  decreasing  downstream  pressure  to 
maintain  mass  conservation  through  the  narrowing  collapsed  section  of  the  tube.  Concurrent  with 
or  following  collapse,  large  or  small  amplitude  oscillations  of  the  tube  walls  may  also  be  seen; 
small  amplitude  wall  oscillations  may  also  occur  in  less  flexible  (uncollapsed)  tubes  transporting 
fluids  at  sufficiently  high  velocities. 

The  oscillatory  wall-fluid  instabilities  that  are  present  in  flexible  tubes  may  be  divided 
into  two  categories;  those  that  require  wall  inertia  (which  we  call  flutter)  and  those  that  do  not. 
The  work  on  oscillatory  instabilities  in  the  absence  of  wall  inertia  is  the  subject  of  a  review  by 
Kamm  and  Pedley  (1989);  we  refer  the  reader  to  their  paper  for  a  complete  treatment  of  and 
reference  list  for  this  subject.  Instabilities  of  this  type  are  commonly  large  amplitude,  lower 
frequency  oscillations  due  to  unsteady  head  loss  resulting  from  the  motion  of  the  separated  flow 
region  downstream  of  a  constriction  in  the  tube.  Flutter  instabilities  are  small  amplitude,  high 
frequency  oscillations  arising  from  the  coupling  of  the  fluid-dynamic  pressure  with  the  compliant 
wall.  Instabilities  similar  in  mechanism  to  Kelvin-Helmholtz  and  water  waves  are  among  these. 
For  the  appearance  of  a  flutter  instability,  a  well  defined  critical  fluid  velocity  must  be  attained. 
Flow  limitation  is  not  a  prerequisite  for  the  instability,  but  flutter  may  be  facilitated  through  the 
decrease  in  tube  cross-sectional  area  and  resulting  increase  in  fluid  flow  speed  that  accompany 
limitation,  and  in  experimental  studies  of  flutter  in  thick-walled  tubes  (Gavriely  et.al.  1989)  and 
wheezing  in  the  lung  (Gavriely  et.al.  1987)  these  oscillatory  phenomena  have  in  fact  been 
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observed  only  in  the  presence  of  flow  limitation. 

§1.3  Flow  limitation  and  flutter  in  the  lung 

Flow  limitation  is  well  documented  in  the  lung,  and  is  described  in  a  review  by  Hyatt 
(1983).  In  healthy  subjects,  flow  limitation  may  be  seen  in  forceful  exhalation,  and  in  lung 
obstructed  patients  the  effort  required  to  cause  flow  limitation  may  be  low  enough  that  the  flow 
rate  in  normal  breathing  may  become  limited.  In  collapsible  tubes  oscillatory  instabilities  are  seen 
following  tube  collapse  and  flow  limitation  (Gavriely  et.al.  1989),  and  it  was  proposed  by 
Grotberg  and  Davis  (1980)  that  wheezing  lung  sounds  may  be  symptomatic  of  such  a  flutter 
instability  in  the  lung  airways.  Experiments  with  forced  expiratory  wheezes  and  lung  preparations 
have  shown  that  flow  limitation  is  a  requisite  for  wheezing  sounds  (Gavriely  et.al.  1987),  that 
these  sounds  are  dependant  on  wall,  rather  than  gas,  properties  (Shabtai  et.al.  1992),  and  that  they 
depend  on  the  flexibility  of  lung  structures  (Gavriely  and  Grotberg  1988).  Thus  a  wall-fluid 
instability  provides  a  likely  explanation  for  the  wheezing  sounds  heard  on  forced  exhalation  and 
in  lung  preparations.  Further,  as  the  spectral  content  and  other  characteristics  of  these  sounds  are 
similar  to  those  of  wheezes  in  patients  with  asthma  and  other  lung  obstructive  diseases  (Gavriely 
et.al.  1987)  and  tube  experiments  (Gavriely  et.al.  1989),  the  sound-producing  mechanism  in  each 
is  likely  to  be  the  same.  In  tube  experiments,  in  particular,  oscillatory  instabilities  may  be 
observed  directly  and  compared  with  theoretical  investigations  of  flutter  (Grotberg  and  Gavriely 
1989)  to  demonstrate  that  this  provides  a  likely  explanation  for  experimental  observations. 

Other  possible  mechanisms  for  the  generation  of  wheezing  lung  sounds  were  considered 
by  Gavriely  et.al.  (1984).  Of  these,  the  only  one  consistent  with  the  experimental  observations 
described  above  (besides  flutter)  is  vortex  resonance,  the  resonance  of  vortices  shed  from  a  sharp 
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edge  or  constriction  with  the  natural  frequency  of  the  compliant  airway  wall,  thus  exciting 
oscillations.  It  appears,  however,  that  there  is  more  evidence  to  support  flutter  as  the  cause  of 
wheezing  lung  sounds. 

§1.4  Organization  of  the  remainder  of  the  dissertation 

In  the  preceding  we  have  outlined  the  experimentally  observed  behavior  of  compliant 
tubes  and  the  evidence  linking  this  to  the  lung  airways.  In  particular,  we  observed  support  for 
the  theory  that  wheezing  lung  sounds  are  symptomatic  of  flutter  of  the  airway  walls.  In  the  next 
chapter  we  review  the  theoretical  work  providing  insight  on  the  problem  of  oscillatory  instabilities 
in  fluid-conveying  compliant  tubes  and  channels.  We  develop  a  model  of  a  collapsed  lung  airway 
or  flexible  tube  in  Chapter  3,  solutions  for  which  are  obtained  in  Chapters  4  and  5.  In  Chapter 
6  we  compare  results  obtained  with  these  models  with  tube  experiments  and  forced  expiratory 
wheezes.  In  Chapter  7  we  conclude  our  analysis  with  a  final  discussion  of  the  results  obtained. 


CHAPTER  2:  LITERATURE  REVIEW 


§2.1  Overview 

Given  the  wide  range  of  physiological  and  engineering  applications  involving  flow 
through  compliant  tubes,  it  comes  as  no  surprise  that  there  is  an  expansive  body  of  literature 
involving  the  theoretical  analysis  of  such  systems.  In  this  chapter  we  review  the  more  relevant 
of  these  works. 

For  flow  through  collapsible  tubes  there  is  a  volume  of  work  with  lumped  parameter  and 
one-dimensional  models.  Lumped  parameter  models  ( e.g .  Conrad  1969;  Katz  et.al.  1969) 
condense  the  global  effects  of  the  tube  geometry  and  fluid  flow  into  an  equation  for  a  simple 
time-dependent  variable,  and  illustrate  the  coupling  of  effects  within  the  system,  but  are  otherwise 
of  limited  application.  One-dimensional  models  consider  the  fluid  variables  and  tube  cross- 
sectional  area  to  be  functions  of  time  and  distance  along  the  tube;  these  are  found  through 
solution  of  the  fluid  equations,  supplemented  by  a  ’tube  law’  relating  the  tube  cross-sectional  area 
to  the  transmural  pressure.  These  models  have  been  applied  to  the  analysis  of  steady  tube 
collapse  (Shapiro  1977)  and  to  unsteady  behavior  (Kamm  and  Shapiro  1979)  and  flow-induced 
oscillations  (Cancelli  and  Pedley  1985)  in  the  absence  of  inertia.  However,  while  these  studies 
have  been  fundamental  in  understanding  the  behavior  of  collapsible  tubes,  they  do  not  include 
investigations  of  the  flutter  instability.  We  therefore  do  not  review  this  field  here,  and  instead 
refer  the  interested  reader  to  the  review  of  Kamm  and  Pedley  (1989). 

While  we  are  specifically  interested  in  studying  the  behavior  of  compliant  fluid- 
transporting  tubes  (or  channels),  there  is  considerable  overlap  between  the  nature  of  the  flutter 
instability  in  such  systems  and  that  for  flows  over  single  compliant  plates.  Our  review  therefore 
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covers  work  on  both  subjects,  and  while  we  use  the  distinction  between  the  geometries  to 
discriminate  between  studies,  the  chapter  sections  are  divided  according  to  other  criteria  that  are 
more  significant  in  determining  the  characteristics  of  the  instabilities  in  the  system.  It  is  perhaps 
not  coincidental  that  these  differences  correspond  roughly  to  the  motivations  for  the  studies. 
Before  continuing  with  the  review  proper  we  briefly  summarize  these  motivations  and  outline  the 
organization  of  the  remainder  of  the  chapter. 

In  §2.2  we  review  some  of  the  earliest  investigations  of  the  stability  of  flows  with 
compliant  boundaries,  which  were  inspired  by  the  experimental  work  of  Kramer  (1960)  on  drag 
reduction  by  use  of  surface  compliance.  The  observation  of  oscillations  in  pipelines  and  other 
structures  spurred  investigation  of  the  stability  of  inviscid  flow  through  tubes  and  over  panels; 
we  consider  this  in  §2.3.  In  §2.4  we  describe  channel  studies  motivated  by  interest  in  both 
engineering  and  physiological  applications,  which  clarify  some  of  the  behavior  seen  previously. 
Some  recent  work  on  the  stability  of  flow  over  a  compliant  plate  and  the  possibility  of  drag 
reduction  by  use  of  a  compliant  surface  are  reviewed  in  §2.5.  Finally  in  §2.6  we  briefly  consider 
some  of  the  other  approaches  in  the  analysis  of  the  behavior  of  systems  involving  flow  past 
compliant  structures,  including  forcing  and  full  numerical  simulations.  In  §2.7  we  conclude  our 
review  in  the  context  of  the  present  work. 

Clearly  the  investigation  of  such  stiff  structures  as  pipelines  will  proceed  from  a  slightly 
different  viewpoint  than  will  work  on  the  more  flexible  tubes  that  are  found  in  physiological 
applications.  Even  for  tubes  that  are  not  for  flow  rates  of  general  interest  susceptible  to  the 
collapse  phenomenon  seen  in  more  flexible  tubes,  however,  the  stability  analysis  of  the  flutter 
instability  may  proceed  in  a  similar  manner  to  that  for  more  flexible  systems,  and  as  a  result 
insight  may  be  gained  from  considering  these  studies. 
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§2.2  Early  work  and  instability  types 

The  pioneering  woric  on  the  interaction  of  flows  with  compliant  structures  is  that  of 
Benjamin  and  Landahl,  who  in  the  early  1960s  considered  the  stability  of  flow  over  a  compliant 
panel.  In  that  the  Tollmien-Schlichting  instability  is  due  to  viscosity,  which  is  confined  to  a 
narrow  boundary  layer  near  the  wall,  Benjamin  (1960)  observed  that  wall  flexibility  should  have 
a  significant  effect  on  the  instability.  His  and  later  work  (Landahl  1962,  Benjamin  1963) 
confirmed  this,  and  also  showed  that  there  are  in  addition  to  the  viscous  fluid  instability  other 
wall-fluid  instabilities  which  arise  as  a  direct  result  of  wall  compliance  and  which  may  exist  in 
the  absence  of  viscosity;  these  include  the  instability  corresponding  to  the  channel  flutter  in  which 
we  are  interested.  Benjamin  and  Landahl  classified  these  instabilities  as  Types  A,  B  or  C  waves 
according  to  the  change  of  the  total  kinetic  and  elastic  energy  in  the  system  necessary  for  their 
growth.  Type  A  instabilities  are  negative  energy  waves  in  the  sense  that  these  energies  decrease 
as  the  amplitude  of  the  wave  increases;  Type  B  are  precisely  the  opposite,  and  for  Type  C  waves 
the  energy  remains  unchanged.  The  stability  of  Types  A  and  B  waves  is  thus  determined  by  the 
net  effect  of  non-conservative  forces  in  the  system  while  Type  C  waves  are  destabilized  by  the 
unidirectional  transfer  of  energy  from  the  fluid  to  the  wall  by  conservative  forces,  and  are 
analogous  to  the  Kelvin-Helmholtz  instability  of  supeiposed  fluids.  Type  C  waves  arise  as  the 
coalescence  of  a  Type  A  and  a  Type  B  wave,  and  are  the  only  type  of  instability  possible  in  a 
conservative  system.  The  Tollmien-Schlichting  instability  may  be  classified  as  Type  A.  As  the 
growth  of  Type  A  waves  requires  a  decrease  in  system  energy  such  waves  will  be  destabilized 
by  wall  damping,  while,  conversely,  Type  B  waves  will  be  stabilized.  Miles  (1957)  demonstrated 
that  for  a  shear  flow  over  a  wavy  wall  there  may  exist  a  component  of  fluid  pressure  in  phase 
with  the  slope  of  the  wall;  this  pressure  component  may  do  work  on,  and  thus  add  energy  to,  the 
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wall.  Through  this  effect,  the  addition  of  viscosity  to  the  problem  may  stabilize  Type  A  waves 
and  destabilize  Type  B  waves.  Dissipative  processes  do  not  affect  the  stability  of  Type  C  waves, 
but  do  cause  a  slight  decrease  in  their  growth  rate. 

Benjamin  and  Landahl  showed  that  as  the  wall  flexibility  is  increased  the  Tollmien 
Schlichting  instability  (TSI)  is  stabilized.  To  postpone  the  onset  of  instability,  however,  it  is  also 
necessary  to  avoid  the  destabilization  of  the  other  wave  types,  so  that  the  compliant  surface  must 
be  chosen  with  care  to  obtain  any  significant  improvement  in  stability  (Landahl  1962).  A  similar 
stabilization  of  the  TSI  was  seen  by  Hains  and  Price  (1962)  for  Poiseuille  flow  through  a  flexible 
channel.  A  good  review  of  this  early  work  appears  in  Benjamin  (1964). 

§2.3  Work  with  finite  length  panels  and  pipes  in  potential  flow 

We  next  consider  work  on  the  stability  of  systems  consisting  of  compliant  panels  or  tubes 
of  finite  length  with  potential  flows.  As  the  flows  are  inviscid  no  information  is  gained  on  the 
behavior  of  the  TSI,  but  the  nature  of  the  wall-fluid  flutter  instability  may  be  examined. 

Weaver  and  Unny  (1970)  considered  the  stability  of  a  system  with  potential  flow  over  a 
two  dimensional  flexible  plate  of  finite  length.  They  found  a  divergence  instability  as  the  flow 
speed  is  increased,  that  is,  an  instability  with  zero  phase  speed;  physically  divergence  may  take 
the  form  of  a  stationary  or  slowly  moving,  growing  wave  train  on  the  flexible  surface.  At  onset 
this  is  a  mode  one  instability,  meaning  that  a  single  wave  spans  the  entire  surface,  but  Weaver 
and  Unny  predicted  that  at  higher  flow  speeds  oscillations  may  be  seen,  followed  by  a  second 
divergence  mode.  However,  in  that  they  used  linear  theory  it  is  not  clear  that  their  conclusions 
regarding  the  behavior  of  the  system  following  the  onset  of  instability  are  valid.  The  addition  of 
damping  to  the  wall  has  little  effect  on  the  instability.  Similar  results  were  seen  by  Komecki 
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(1974),  Komecki  et.  al.  (1976)  and  Ellen  (1973). 

The  results  for  finite  length  tubes  conveying  potential  flows  are  similar  to  those  for 
panels.  Weaver  and  Unny  (1973)  and  Komecki  (1974)  examined  the  behavior  of  such  systems 
and  determined  that  the  instability  is  again  divergent  at  criticality,  with  flutter  occurring  only  at 
higher  flow  speeds,  if  at  all  (though,  again,  the  use  of  linear  theory  calls  into  question  conclusions 
drawn  following  the  onset  of  instability).  This  flutter  is  a  ’coupled  mode’  instability  arising  from 
the  coalescence  of  two  divergence  modes;  the  manner  in  which  this  coupling  takes  place  was 
elucidated  by  Paidoussis  and  Issid  (1974),  who  considered  the  problem  with  a  more  accurate  pipe 
model.  Depending  on  the  choice  of  wall  characteristics  the  flutter  instability  may  or  may  not  be 
present,  though  when  present  it  is  always  preceded  by  divergence.  Paidoussis  and  Issid ’s  work 
also  includes  a  review  of  the  literature  on  the  problem  of  the  stability  of  flow  in  flexible  pipes. 

§2.4  Work  with  channel  flow 

A  channel  may  be  seen  as  a  two  dimensional  model  of  a  tube,  but  there  are  additional 
physical  considerations  that  suggest  that  the  channel  geometry  is  worthy  of  attention.  Weaver  and 
Paidoussis  (1977)  observed  that  the  flutter  instability  seen  in  tubes  is  not  the  classical  shell  flutter 
mode  in  which  the  tube  flattens  alternately  in  the  two  perpendicular  directions  normal  to  the  tube 
axis,  but  rather  a  ’flapping  flutter’  that  occurs  after  the  tube  has  flattened  somewhat  to  become 
oval.  This  flapping  flutter  then  involves  the  oscillation  of  the  longer  opposing  walls  of  this  oval 
either  in  or  out  of  phase  with  one  another,  and  its  nature  suggests  that  it  might  be  appropriate  to 
model  the  tube  as  a  channel.  A  similar  conclusion  may  be  reached  on  consideration  of  the  nature 
of  the  collapse  and  oscillation  of  the  flexible  tubes  appearing  in  physiology;  these  tubes  are  more 
flexible  than  those  considered  by  the  authors  mentioned  in  §2.3,  so  that  there  may  be  dramatic 


10 


changes  in  their  cross-sectional  geometry,  as  noted  in  §1.2.  In  this  case  a  circular  tube  model  is 
not  representative  of  the  actual  tube  geometry,  and  it  may  be  more  appropriate  to  instead  model 
the  central  (flat)  section  of  the  collapsed  tube  (which  is  that  which  oscillates)  as  a  channel. 

Motivated  by  their  observation,  Weaver  and  Paidoussis  considered  a  channel  with  flexible 
walls  conveying  an  inviscid  flow.  They  used  two  wall  models,  one  with  inertia,  bending  stiffness 
and  damping,  and  another  consisting  of  two  parallel  plates  with  constrained  edges  but  infinite 
length  (so  that  the  channel  was  in  this  case  essentially  three  dimensional,  but  without  side  walls). 
In  either  case  they  found  a  divergence  instability  for  both  the  symmetric  and  antisymmetric  flutter 
modes  (in  which  the  channel  walls  oscillate  out  of  or  in  phase  with  one  another,  respectively), 
followed  by  coupled  mode  flutter.  Again,  only  linear  theory  was  considered,  so  that  their 
demonstration  of  flutter  remains  inconclusive.  They  compared  these  results  with  experiments  and 
found  qualitative  agreement.  (It  should  be  noted  that  Weaver  and  Paidoussis  used  the  method  of 
images  to  find  the  fluid  velocity;  however,  while  this  results  in  an  expression  involving  an  infinite 
sum,  they  included  only  the  first  two  terms,  so  that  their  solution  allows  cross -flow  at  the  walls 
of  the  channel.)  Matsuzaki  and  Fung  (1977)  considered  a  similar  problem,  with  inviscid  flow 
through  a  channel  with  a  flexible  section  of  finite  length,  modelled  using  the  von  K&rm&n  plate 
equations,  and  also  found  a  divergence  instability. 

Grotberg  and  Davis  (1980)  considered  the  channel  flow  problem  (with  flexible  walls  of 
infinite  length)  as  a  model  of  a  collapsed  lung  airway,  theorizing  that  wheezing  lung  sounds  may 
be  symptomatic  of  a  flutter  oscillation  in  the  airway  walls.  They  used  two  wall  models  that 
included  inertia  and  either  bending  stiffness  or  elastance,  respectively,  and  found  a  flutter 
instability  unless  wall  damping  was  also  included  in  the  wall,  in  which  case  instability  appeared 
at  lower  flow  speeds  and  was  divergent.  This  is  consistent  with  the  results  of  Weaver  and 
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Pai'doussis  (1977),  as  the  latter  considered  only  damped  walls.  Grotberg  and  Davis  showed  also 
that  results  obtained  in  the  two  dimensional  geometry  are  consistent  with  one  dimensional  studies 
(e.g.  Shapiro,  1977).  The  drop  in  critical  flow  speed  with  the  addition  of  damping  was  first 
explained  by  Landahl  (1962),  who  showed  that  without  wall  damping  there  is  preceding  the 
appearance  of  instability  a  state  of  unstable  equilibrium  (neutral  stability)  that  is  broken  through 
the  addition  of  damping,  thereby  reducing  the  critical  flow  speed.  A  similar  drop  in  critical  flow 
speed  may  be  obtained  by  including  fluid  viscosity  in  the  absence  of  wall  damping;  this  results 
in  the  appearance  of  a  flutter  instability.  This  observation  was  made  for  flow  over  a  compliant 
plate  by  Carpenter  (1984),  and  may  be  seen  for  channel  flow  with  the  Grotberg  and  Reiss  (1984) 
model  (c./.  below). 

The  appearance  of  divergence  only  in  systems  with  wall  damping  appears  at  first  sight 
to  contradict  the  results  seen  in  §2.3,  in  which  irrespective  of  whether  the  walls  were  damped  the 
instability  was  divergent  at  onset.  This  difference  is,  however,  due  to  another  fundamental 
difference  between  the  systems  being  considered  in  §2.3  and  above,  namely  the  finite  length  of 
the  flexible  section.  Lucey  and  Carpenter  (1992,  1993)  addressed  this  issue  when  considering 
flow  over  a  finite  length  compliant  plate,  and  showed  that  in  the  case  of  a  finite  geometry  there 
are  always  end  effects  that  influence  the  system  in  the  same  manner  as  damping  in  the  infinite- 
length  case. 

In  all  of  the  analyses  thus  far  considered,  a  divergence  instability  is  seen  at  criticality. 
However,  in  physical  systems  flutter  is  observed;  to  address  this  discrepancy,  Grotberg  and  Reiss 
(1982,  1984)  added  a  hydraulic  friction  term  to  the  potential  flow  through  a  channel  to 
approximate  the  effect  of  viscosity.  (Grotberg  and  Shee  (1985)  showed  that  this  is  consistent  with 
the  full  Orr-Sommerfeld  system  when  viscosity  and  wall  damping  are  taken  to  zero  and  their  ratio 


12 


is  0(1).)  They  found  that  the  addition  of  fluid  friction  results  in  the  reappearance  of  the  flutter 
instability,  and  through  a  weakly  non-linear  analysis  showed  also  that  the  bifurcation  to  flutter 
is  supercritical  and  hence  stable.  The  recovery  of  a  flutter  instability  with  the  addition  of 
viscosity  is  due  to  the  fact  that  viscosity  stabilizes  Type  A  wall-fluid  instabilities  (as  noted  in 
§2.2)  while  destabilizing  those  of  Type  B.  The  introduction  of  wall  damping  to  the  inviscid 
system  destabilizes  the  Type  A  wave,  leading  to  a  divergence  instability;  the  addition  of  viscosity 
may  restabilize  this  wave  (and  destabilize  the  Type  B  wave),  resulting  in  an  oscillatory  instability. 
The  Grotberg  and  Reiss  channel  model  was  reconsidered  with  the  addition  of  fluid  compressibility 
by  Grotberg  and  Shee  (1985),  who  found  generally  slight  corrections  to  the  incompressible  case. 
To  show  that  the  wheezing  lung  sounds  that  Grotberg  and  Davis  (1980)  first  set  out  to  describe 
are  in  fact  due  to  flutter  of  the  airway  walls,  Grotberg  and  Gavriely  (1989)  compared  the  results 
obtained  with  the  Grotberg  and  Reiss  model  with  experiments  involving  flow  through  flexible 
tubes,  and  found  good  agreement  between  theory  and  experiments. 

Webster  et.  al.  (1985)  experimentally  and  theoretically  examined  flow  through  a  model 
of  the  trachea,  a  three  sided  duct  covered  with  a  compliant  tensioned  membrane.  Their  theoretical 
model  was  a  two-dimensional  channel  with  one  compliant  boundary;  the  flow  was  taken  to  be 
potential  and  the  wall  undamped.  Due  to  the  singular  influence  of  damping  on  this  type  of 
system,  the  relevance  of  this  analysis  to  physical  situations  might  be  questioned,  but  they  reported 
good  agreement  with  their  experiments. 

§2.5  Additional  work  on  flow  over  compliant  plates 

The  possibility  of  the  stabilization  of  the  Blasius  boundary  layer  by  a  compliant  plate  was 
reconsidered  by  Carpenter  and  Garrad  (Garrad  and  Carpenter  1982,  Carpenter  and  Garrad  1985, 
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1986).  They  considered  flow  over  a  flexible  plate  with  linear  wall  properties,  and  allowed  for 
the  effect  of  a  ’fluid  substrate’  in  a  reservoir  behind  the  plate.  Garrad  and  Carpenter  (1982) 
examined  a  compliant  plate  of  finite  length,  with  an  inviscid  substrate  fluid  of  infinite  depth,  and 
a  potential  flow  with  corrections  to  simulate  the  effect  of  a  laminar  or  turbulent  boundary  layer. 
Using  the  Galerkin  method,  they  found  for  the  laminar  case  a  divergence  instability,  and  for  the 
turbulent  case  flutter.  The  persistence  of  the  divergence  instability  here  even  with  the  inclusion 
of  viscosity  is  reasonable,  as  the  flexible  surface  had  finite  length,  which  may  promote  the 
divergence  instability. 

Carpenter  and  Garrad  (1985,  1986)  subsequently  performed  an  in  depth  analysis  of  the 
linear  stability  of  the  infinite  length  system  both  with  and  without  a  viscous  substrate  fluid, 
considering  both  the  Tollmien-Schlichting  and  wall-fluid  instabilities  (in  the  1985  and  1986 
papers,  respectively).  Approximate  analytical  methods  were  used  to  examine  how  different 
physical  effects  alter  the  stability  of  the  different  waves,  thereby  illustrating  the  behavior  of  the 
wave  types  of  Benjamin  and  Landahl.  The  linear  stability  analysis  of  the  TSI  was  carried  out 
numerically  to  account  for  the  full  effects  of  viscosity,  using  a  shooting  algorithm.  For  the  flutter 
instability  this  numerical  solution  was  supplemented  by  an  analytical  analysis  using  Benjamin’s 
(1963)  approximation  to  the  pressure  fluctuation  caused  by  a  boundary  layer.  They  then 
addressed  the  question  of  whether  Kramer’s  (1960)  experiments  could  have  shown  significant 
postponement  of  the  transition  to  turbulence,  and  concluded  that  this  was  in  fact  a  possibility. 
Their  1985  paper  also  includes  a  review  of  the  literature  on  the  theoretical  and  experimental 
stability  analyses  for  flow  over  a  flat  plate.  Other  such  reviews  are  those  of  Gad  el  Hak  (1986) 
and  Riley  et.al.  (1988). 

More  recently,  Caipenter  and  Gajjar  (1990)  formulated  a  general  multi-deck  asymptotic 
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theory  for  the  stability  of  the  Type  B  (travelling  wave)  flutter  instability  in  flow  over  a  compliant 
plate  of  infinite  length.  This  is  significantly  less  computationally  intensive  than  the  solution  of 
the  full  Orr-Sommerfeld  system,  but  is  limited  to  this  flutter  instability  and  the  case  in  which  the 
critical  layer  (where  the  phase  speed  of  the  disturbance  becomes  equal  to  the  flow  speed  of  the 
base  flow)  is  distinct  from  the  viscous  boundary  layer  at  the  wall. 

Additional  investigations  of  the  linear  stability  of  a  flow  over  a  compliant  surface  are 
those  of  Duncan  et.al.  (1985)  and  Evrensel  and  Kalnins  (1985,  1988).  Duncan  et.al.  considered 
a  potential  flow  with  a  phase  correction  to  the  perturbation  pressure  to  simulate  a  boundary  layer, 
while  Evrensel  and  Kalnins  in  1985  considered  a  potential  flow,  and  in  1988  a  fully  viscous 
profile.  The  (linear)  stability  calculations  in  each  case  differ  from  those  of  Carpenter  and  Garrad 
(1985,  1986)  primarily  in  the  wall  model  chosen;  these  authors  consider  the  wall  to  be  an  elastic 
or  viscoelastic  solid  of  finite  depth  and  solve  the  appropriate  equations  for  the  behavior  of  such 
a  material  to  determine  the  wall  motion.  Pierucci  and  Morales  (1990)  considered  the  equivalent 
problem  for  the  TSI  in  plane  Poiseuille  flow  with  elastic  (undamped)  walls,  but  their  conclusions 
are  questionable  due  to  their  assumption  that  the  critical  wave  number  remains  unchanged  with 
the  introduction  of  wall  compliance. 

§2.6  Other  approaches 

In  addition  to  the  stability  analyses  described  in  the  preceding  sections,  other  methods 
have  been  used  to  obtain  a  better  understanding  of  the  nature  of  the  coupling  of  flow  to  compliant 
structures  and  relevant  applications.  In  this  section  we  briefly  review  some  of  these  studies. 

Walsh  et.al.  (1991)  modelled  expiratory  flow  in  the  trachea  with  a  geometry  consisting 
of  two  concentric  cylinders,  the  inner  of  which  has  an  undamped,  untethered  flexible  wall.  Using 
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curved-shell  flutter  theory  with  one-dimensional  inviscid  fluid  equations,  they  considered  the  axial 
displacements  of  the  wall  to  be  much  larger  than  the  lateral  displacements.  With  these 
assumptions,  they  found  a  long-wave  flutter  instability. 

Brazier-Smith  and  Scott  (1984)  considered  the  linear  analysis  of  an  inviscid  flow  over  a 
temporally  forced  undamped  plate  with  bending  stiffness.  As  the  system  is  both  undamped  and 
forced  it  is  unstable,  and  the  resulting  instability  was  studied  to  determine  whether  it  is  absolute 
or  convective.  Carpenter  and  Garrad  (1986)  compared  their  results  with  those  of  Brazier-Smith 
and  Scott  by  calculating  the  group  velocity  of  the  disturbance  to  the  system  as  the  derivative  of 
the  phase  speed  with  respect  to  wave  number,  and  found  good  agreement.  Lucey  and  Carpenter 
(1992)  presented  a  similar  but  numerical  analysis  for  an  inviscid  flow  over  a  forced  compliant 
plate  of  finite  length.  They  were  concerned  exclusively  with  the  divergence  instability,  which  was 
isolated  through  consideration  of  the  energy  of  the  system. 

A  full  numerical  solution  of  the  Navier  Stokes  and  compliant  wall  equations  was  obtained 
by  Domaradzki  and  Metcalfe  (1987)  for  flow  over  a  flexible  plate,  and  used  to  determine  the 
nauire  of  the  distribution  and  dissipation  of  energy  for  the  different  wave  types  of  Benjamin  and 
Landahl.  For  Type  A  waves  the  energy  dissipation  has  a  local  maximum  near  the  wall,  while  for 
Type  B  waves  there  is  little  energy  activity  there.  Type  A  waves  are  stabilized  primarily  by  the 
boundary  layer  region  of  fluid,  while  wall  dissipation  is  the  main  stabilizing  factor  for  Type  B 
waves.  These  observations  are  clearly  representative  of  the  general  characterizations  of  the  given 
wave  types. 

Sen  and  Arora  (1988)  approached  the  problem  of  flow  over  a  compliant  wall  by  using 
a  kinematic  formulation  of  the  wall  boundary  conditions  that  permits  analysis  of  the  flow  without 
specification  of  the  specific  material  properties  of  the  wall.  These  properties  may  then  be  back 
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calculated  given  the  wall  dynamics.  By  thus  separating  the  nature  of  the  wall  model  from  the 
actual  analysis  of  the  flow,  they  proposed  that  walls  with  optimum  stabilization  capabilities  could 
be  determined.  A  general  understanding  of  the  instability  for  flows  with  compliant  boundaries 
was  also  the  objective  of  Yeo  and  Dowling  (1987),  who  reconsidered  many  of  the  fundamental 
stability  theorems  for  inviscid  flow  over  a  rigid  surface  to  extend  their  validity  to  the  case  of  a 
compliant  wall. 

Recently,  there  have  been  a  few  non-linear  studies  of  the  TSI  for  Poiseuille  flow  in  a 
channel.  Rotenberry  and  Saffman  (1990)  undertook  a  weakly  non-linear  analysis  of  the  instability 
and  derived  a  Ginzburg-Landau  equation  for  the  amplitude  of  the  bifurcating  solution.  The  walls 
considered  were  simple  elastic  membranes  (with  allowance  for  damping),  and  they  found  that 
when  these  are  sufficiently  flexible  the  bifurcation  to  the  TSI  goes  from  being  subcritical,  as  is 
the  case  for  a  channel  with  rigid  walls,  to  supercritical.  Rotenberry  and  Saffman  also  proved  an 
extension  of  Squire’s  theorem  for  Poiseuille  flow  in  a  compliant  channel.  Rotenberry  (1992) 
considered  the  same  channel  system  with  the  addition  of  tension  to  the  channel  walls,  and  inferred 
the  criticality  of  the  bifurcation  from  the  shape  of  the  curve  of  disturbance  energy  as  a  function 
of  Reynolds  number.  He  concluded  that  while  sufficiently  flexible  walls  may  render  the 
bifurcation  supercritical,  the  magnitude  of  the  supercritical  branch  is  in  fact  very  small,  so  that 
for  finite  but  small  distances  from  the  base  Poiseuille  solution  the  bifurcation  appears  subcritical. 
Thus  wall  flexibility  does  not  quantitatively  change  the  character  of  the  TSI. 

§2.7  Conclusions 

These  studies  show  that  the  effects  of  wall  and  fluid  damping  are  significant  in 
determining  the  nature  of  the  flutter  instability,  and  that  their  influence  is  often  singular  in  nature. 
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Thus  our  model  is  developed  to  accurately  include  these  effects,  and  is  two-dimensional  to  allow 
for  the  effect  of  transverse  variation  in  the  flow  speed  in  the  channel  (which  has  not  been 
previously  considered  in  channel  studies  of  the  flutter  instability).  The  similarity  between 
previous  results  for  and  methods  used  in  the  study  of  more  rigid  (uncollapsed)  pipes  and 
collapsible  tubes  suggests  that  the  channel  model  we  develop  will  also  be  applicable  to  a  wide 
range  of  systems,  though  we  are  specifically  interested  in  its  application  to  the  lung  airways  and 
wheezing  lung  sounds. 


CHAPTER  3:  MODEL  DERIVATION 


§3.1  Derivation  of  nondimensional  equations  and  boundary  conditions 

To  model  a  partially  collapsed  flexible  tube  we  consider  a  two-dimensional  channel  with 
compliant  walls;  through  this  flows  an  incompressible  Newtonian  fluid.  In  this  chapter  we  derive 
the  linear  stability  equations  for  this  system  with  a  basic  state  consisting  of  a  flat  walled  channel 
and  a  given  laminar  flow. 
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Figure  3.1  Flexible  walled  channel  system.  Walls  are  spring  supported  and 
damped,  and  a  coordinate  system  is  fixed  to  the  channel  midline; 
dimensions  are  as  shown. 

The  system  under  consideration  is  shown  in  figure  3.1.  The  walls  of  the  channel  are 
damped  flexible  plates  supported  by  an  array  of  springs  attached  to  a  rigid  backing;  it  is  assumed 
that  this  continuous  array  of  linear  springs  models  the  elastance  of  the  tube  and  the  material 
supporting  it.  We  define  a  coordinate  system  centered  in  the  channel  with  axial  and  transverse 


18 


19 


coordinates  x*  and  z*.  As  indicated,  the  half  channel  width  is  b‘,  so  that  the  undisturbed  location 
of  the  walls  is  z*  =  ±  b“.  Denoting  the  (x\z* )  location  of  a  point  on  the  lower  wall  as  U*  = 
(U*,W* ),  the  nonlinear  von  Kdrmdn  plate  equations  governing  the  normal  and  tangential  motion 
of  the  wall  are 


,,  \2D' (dU'  d2W'  d2U'dW'  ,dW\id2W\ 

M  ipr  +  G  -Sr  -  — (1F S*  *  lF~aF  *  ‘iF*  TS^ 

-  T^W*  -  E' (W*  +  b')  +  (S^-S^r)  =  0 
dx'*  3x*2  (3.1) 


and 


p + 

Kw  ar2 

yr  aw* 


3r*  h *2  9x’ 

.)  +  (S'T  -  $£■)  =  0 


(3.2) 


3x*2  3x* 

(Sapir  and  Reiss  1979),  where  pw*  is  the  density  of  the  wall  material,  h*  the  wall  thickness,  G* 
and  H*  the  coefficients  of  normal  and  tangential  damping  in  the  wall,  D*  the  flexural  rigidity  of 
the  plate,  T*  the  imposed  longitudinal  tension,  E*  the  spring  constant  (elastance)  of  the  spring 
supports,  and  S*N  and  S*T  the  normal  and  tangential  stresses  on  the  wall,  respectively,  with 
subscript  ’ext’  indicating  stresses  acting  from  the  region  external  to  the  wall.  As  written,  these 
equations  apply  to  the  lower  wall  of  the  channel.  In  the  development  that  follows,  we  for 
simplicity  treat  only  the  equations  for  the  lower  wall;  in  all  cases  similar  equations  apply  at  the 
upper  wall.  In  equations  (3.1,  2)  and  those  following,  a  superscripted  asterisk  indicates  a 
dimensional  variable. 


The  motion  of  the  fluid  is  given  by  the  Navier-Stokes  equations,  in  vector  form 
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p»  {  d“*  +  u’-Vu*)  =  -V*p*  +  v*V,2«* 

dr 


(3.3) 


and 


(3.4) 


V*>«*  =  0 

where  u*  =  (u\w‘)  is  the  fluid  velocity  vector  in  the  x*-z*  plane,  V*  the  two-dimensional  (2D) 
gradient  operator,  V*2  the  2D  laplacian,  and  p*.  p*  and  v*  the  fluid  density,  pressure  and  viscosity, 
respectively.  The  fluid  must  satisfy  the  kinematic  and  no  slip  conditions  at  the  wall, 


n-u *  =  n 


dU' 


and  x  ■«*  =  x 


(3.5) 


dr  dr 

where  n  and  x  are  unit  vectors  normal  and  tangential  to  the  wall.  Hereafter  we  denote  partial 
differentiation  with  subscripts. 

Equations  (3.1-5)  are  nondimensionalized  using  a  velocity  scale  fi,  length  scale  L,  time 
scale  T0  =  L  /  fl  and  pressure  scale  k  =  p*  fl2;  we  generally  take  L  =  b*  (the  half  channel  width) 
and  fl  =  (E*  b*  /  p)1/2,  an  elastic  wave  speed  of  the  wall.  Introducing  the  nondimensional  variables 
U  =  (U,W)  =  U*  /  L,  u  =  (u,w)  =  u*  /  fl,  t  =  t*  /  T0,  (x,  z)  =  (x‘,  z  )  /  L  and  p  =  p’  /  it  (and 
noting  that  the  nondimensional  stresses  are  SN,T  =  S‘N,T  /  jc),  equations  (3.1)  and  (3.2)  become 


MW„  +  2 GW,  -  d(U MW„  +  UaWx  +  (wg2w„) 

+  -  rw  +  £(W  +  ! 1)  +  (SN-S£cr)  =  0 

JBUOC  XX  ¥ 


(3.6) 


and 


MUn  +  2HUt  -  d(Ua  +  W^W;)  +  (Sr  -  5^)  =  0.  (3-7> 

Nondimensional  parameters  appearing  in  (3.6)  and  (3.7)  are  M  =  (pw*  h*  /  p*  L),  G  =  M  L  G* 
/  fi,  H  =  M  L  H*  /  fi,  d  =  12  D*  /  h*2  L  p*  Q2,  B  =  D*  /  L3  p*  fl2,  T  =  T*  /  L  p*  Q2,  and  E  = 
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E*  L  /  p*  02  (  =  1  for  the  scaling  indicated  above).  These  parameters  are  respectively  the  wall  to 
fluid  mass  ratio,  and  the  non-dimensional  wall  characteristics:  normal  and  tangential  damping, 
bending  stiffness,  flexural  rigidity,  imposed  longitudinal  tension,  and  elastance. 

The  stresses  in  (3.6)  and  (3.7)  are  found  from  consideration  of  the  fluid  stress  tensor  o, 
in  terms  of  which  SN  =  n  •  o  •  n  and  ST  =  x  •  o  •  n.  Assuming  deflections  of  the  plate  to  be 
small,  these  are 


S"  =  ~P  * 


and  ST  = 


T  _  +  Wx) 


R 


(3.8) 


where  ^  =  (1  L  /  v  is  the  Reynolds  number.  We  have  chosen  the  notation  Rw  to  emphasize  that 
this  Reynolds  number,  due  to  our  choice  of  fl  as  a  wall  elastic  wave  speed,  is  related  to  the  wall 
elastance. 


It  is  convenient  to  rewrite  the  Navier-Stokes  equations  (3.3)  and  (3.4)  in  stream  function 
form  to  eliminate  the  fluid  pressure.  Nondimensionalizing  as  above  and  introducing  a  stream 
function  'F,  where  (u,  w)  =  (4**,  -4'*)  (so  that  the  continuity  equation  (3.4)  is  automatically 
satisfied),  equation  (3.3)  becomes  the  voiticity  equation 


V24>  -  _La24*  =  4*  V24>  -  4#  V24*  , 
'  R  x 


(3.9) 


where  A2  is  the  two-dimensional  biharmonic  operator,  (d*  +  2  9X2  9Z2  +  d*).  The  unit  vectors  n 
and  T  in  (3.5)  are 


n  - 


1 


(.K  +  l) 


1/2 


(-W  ,  1)  and 


t  = 


(w;  +  \)m 


.(i.  K ), 


(3.10) 


so  that  on  introduction  of  the  stream  function,  the  kinematic  and  no  slip  conditions  (3.5)  become 
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-WaHf,-V,  --WM  +  W'  and  (3.ii) 

¥  -  Wx'¥x  =  Ut  +  Wt Wt ,  at  z  =  ±  1 . 

§3.2  System  linearization  and  simplification 

A  basic  solution  for  the  system  defined  by  equations  (3.6-7),  (3.9)  and  (3.11)  has  flat 
walls  (U  =  0,  W  =  ±  1)  and  a  given  stream  function  ¥0  corresponding  to  the  flow  profile  being 

considered  (discussed  below).  Implicit  in  this  base  state  is  the  assumption  of  a  pressure  gradient 

external  to  the  flexible  walls  to  maintain  the  flat  walled  configuration.  We  assume  that  the  base 
flow  is  parallel,  that  is,  that  %  is  a  function  of  z  only,  and  write  d¥0  /  dz  =  u0(z),  where  u^z) 
is  the  base  axial  flow  velocity  in  the  channel. 

To  determine  the  stability  of  this  base  state  we  use  linear  stability  theory  and  introduce 
perturbations  (U'.W')  and  \|/  to  the  base  wall  position  and  stream  function  respectively.  In  the 
usual  manner,  we  write  these  disturbance  quantities  in  normal  modes,  so  that 

(U,W)  *  (0,-1)  +  (U',W')  =  ( 0,  - 1 )  +  (x,tD)exp(i*(x  -  ct))  (3  12) 

¥  =  ¥0  +  y'  =  y0  +  (p(z)exp(ifc(x  -  ct). 

Plugging  these  into  the  governing  equations  and  linearizing  in  perturbation  quantities,  linearized 
stability  equations  are  obtained.  The  vorticity  equation  (3.9)  becomes  the  Orr-Sommerfeld 
equation, 

k(u0  -  c)(<p"  -  k2y)  -  ku0"i p  +  -  2&2cp"  +  &4<p)  =  0,  (3.13) 

Rw 

and  the  kinematic  and  no  slip  conditions  (3.10)  become 

-ikc%  =  <p'  +  u0' co  and  <p  =  (c  -  u0)co,  (3.14) 

where  Uo  and  <p  are,  after  linearization,  evaluated  at  z  =  ±  1.  Using  (3.14)  to  eliminate  the  wall 
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variables  from  (3.6)  and  (3.7),  we  obtain  after  substituting  (3.8)  for  the  fluid  stresses  (and 
eliminating  the  pressure  in  (3.8)  using  (3.3)) 

(-( kc)2M  -  likcG  +  Bk 4  +  Tk2  +  E)kq>  + 

k(c  -  u0)(u0'<. p  +  (c  -  u0)<p')  -  jj — —OP'"  ~  3*V)  =  0  *  * 

Kw 

and 


(~(kc)2M  -  2ikcH  +  dk2)((u0  -  c)<p'  -  <pu0')  - 

ikc(u0  -  c) 


(9"  +  ^9)  =  0. 


(3.16) 


Again,  consistent  with  linearization,  (3.15)  and  (3.16)  are  evaluated  at  z  =  -1. 

By  assuming  disturbances  to  be  divisible  into  ’symmetric’  and  ’antisymmetric’  modes, 
we  may  restrict  attention  to  the  lower  half  of  the  channel.  For  symmetric  disturbances  the 
channel  walls  oscillate  out  of  phase  with  one  another  (in  a  varicose  shape)  so  that  the  disturbance 
stream  function  must  satisfy  conditions  of  no  cross  flow  at  the  midline;  for  antisymmetric 
disturbances  the  walls  are  in  phase  (a  sinuous  shape)  and  we  require  no  axial  flow  at  the  midline. 
These  conditions  are 


9  =  9"  =  0  (symmetric)  or  9'  =  9'"  =  0  (antisymmetric).  (3.17) 

Equation  (3. 13)  with  boundary  conditions  (3.15-17)  constitutes  an  eigenvalue  problem  for 
the  disturbance  stream  function  9  and  eigenvalue  c.  The  real  and  imaginary  parts  of  c  are  the 
phase  speed  and  growth  rate,  respectively,  of  the  disturbance;  instability  thus  sets  in  when  the 
imaginary  part  of  c  becomes  positive,  and  it  is  of  interest  to  determine  the  conditions  for  which 
this  occurs.  The  sign  of  the  real  part  of  c  determines  whether  the  disturbance  travels  upstream 
(negative  sign)  or  downstream  (positive  sign),  and  the  product  of  the  wave  number  k  and  the  real 
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part  of  c  (over  2k)  gives  the  frequency  of  the  oscillations.  For  simplicity  we  assume  that 
tangential  wall  inertia  and  damping  may  be  neglected,  that  is,  that  we  may  discard  the  first  two 
terms  in  the  tangential  boundary  condition  (3.16);  comparison  of  the  numerical  solution  with  the 
resulting  modified  boundary  condition  with  those  using  the  full  condition  shows  the  two  to  be 
indistinguishable.  If  the  tangential  stress  is  also  omitted  from  equation  (3.16)  we  recover  the  no 
slip  condition  used  by  Benjamin  (1960)  and  others,  which  does  not  include  the  effect  of 
horizontal  wall  motion. 

§  3.3  Specification  of  the  base  flow  profile 

To  complete  the  specification  of  the  Orr-Sommerfeld  system  derived  in  §3.2,  a  base  flow 
profile  must  be  chosen.  As  a  solution  for  flow  in  a  channel,  this  will  be  either  a  developing 
(boundary  layer)  or  fully-developed  (Poiseuille)  flow;  we  are  interested  in  modelling  short  tubes, 
and  so  are  concerned  with  the  developing  flow  profile.  However,  there  is  no  closed  form  solution 
for  the  corresponding  stream  function;  as  a  result  there  are  a  large  number  of  approximate 
solutions.  These  may  be  divided  into  four  general  categories:  integral  methods  (Schiller  1922; 
Campbell  and  Slattery  1963);  models  linearizing  the  inertial  terms  in  the  Navier-Stokes  equations 
(Langhaar  1942;  Sparrow  et.al.  1964);  numerical  solutions  of  the  boundary  layer  (Bodoia  and 
Osterle  1961)  or  full  equations  (Brandt  and  Gillis  1964;  Wang  and  Longwell  1964);  and  axially 
patched  or  asymptotic  solutions  (Schlichting  1934;  Van  Dyke  1970;  Wilson  1970, 1971).  A  good 
summary  of  these  categories  and  solution  methods  appears  in  Schmidt  and  Zeldin  (1969). 

We  approximate  the  developing  profile  in  two  ways:  first,  by  using  a  uniform  (plug) 
profile,  and  second,  by  solving  the  boundary  layer  equations  using  either  finite  difference  or 
perturbation  methods.  For  the  flow  regimes  in  which  we  are  interested  the  flow  is  uniform  over 
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most  of  the  channel,  with  nonuniformity  being  restricted  to  a  small  region  (boundary  layer)  near 
the  wall.  The  flow  is  thus  quantitatively  similar  to  a  plug  flow  over  most  of  the  channel  width, 
so  that  it  may  be  hoped  that  meaningful  predictions  may  be  obtained  from  the  analysis  of  this 
much  simpler  flow;  this  analysis  appears  in  Chapter  4.  Of  course,  the  results  so  obtained  must 
be  evaluated  carefully,  as  the  plug  flow  does  not  satisfy  the  no  slip  conditions  at  the  wall. 

Our  finite  difference  solution  uses  a  second  order  scheme  with  a  nonunifoim  grid  to  solve 
the  boundary  layer  equations;  full  details  of  the  method  appear  in  Appendix  A.  This  solution  is 
able  to  provide  the  flow  profile  at  any  location  in  the  channel  up  to  and  including  the  point  of 
fully  developed  flow,  but  more  difficult  to  use  in  the  Orr-Sommerfeld  system  as  any  change  in 
the  velocity  of  the  base  flow  requires  that  the  flow  profile  be  recalculated.  However,  the 
numerically  determined  profile  is  useful  for  evaluating  the  perturbation  solution,  which  we  use 
in  our  stability  calculations. 

We  follow  Schlichting  (1934)  in  the  development  of  a  perturbation  solution  for 
developing  channel  flow.  A  uniform  flow  is  assumed  at  the  channel  inlet  that,  as  it  progresses 
down  the  channel,  assumes  a  profile  with  growing  boundary  layers  at  either  wall.  For  short 
distances  from  the  channel  inlet  these  boundary  layers  form  in  the  same  manner  as  the  Blasius 
solution  for  flow  over  a  flat  plate,  with  the  caveat  that  mass  conservation  requires  that  the  flow 
in  the  central  region  of  the  channel  accelerate  to  compensate  for  the  decreased  flow  at  the  walls. 
Thus  for  the  upstream  portion  of  the  channel  the  difference  between  the  channel  and  flat  plate 
geometries  may  be  incorporated  as  a  perturbation  to  the  Blasius  profile. 

We  therefore  model  the  channel  flow  using  the  steady  boundary  layer  equations,  in 


stream-function  form 
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-  vcv,v  -  u-v;.  ♦  vv;,,..  <318) 

Here  \|/‘  is  the  dimensional  form  of  the  base  flow  stream  function  ¥<,,  0*  (x*)  is  the  velocity  field 
outside  of  the  boundary  layer  and  a  new  transverse  coordinate  y*  =  b*  +  z*  has  been  introduced 
to  mimic  the  transverse  coordinate  usually  used  for  flow  over  a  single  plate.  The  stream  function 
y*  is  written  as  a  perturbation  expansion  that  will  at  leading  order  recover  the  Blasius  solution, 


y*  -  S'b'Uf^r |)  +  e2/^)  +  e3/3(ri)  +  •••).  ^3-19' 

where  e  s  (v*  x*  /  b*2  S*  )1/2,  S*  is  the  velocity  of  the  initial  (uniform)  flow,  and  t|  =  y*  /  (e  b* ) 
is  the  similarity  variable  for  flow  over  a  flat  plate.  The  free  stream  velocity  U  (x)  may  be 
determined  from  mass  conservation,  as  follows.  For  the  flow  rate  at  the  inlet  to  be  equal  to  that 
at  any  point  downstream,  we  must  have  S*  b*  =  U*(b*  -  S/),  where  81*  is  the  displacement 
thickness  of  the  boundary  layer.  By  definition 


b • 

5;  U*  *  J(t7*  -  u')dy'  =  zb'  J  (t/*  -  $•(/,'  +  e/2'  +  \,  (3-2°) 

0  0 

where  r),  is  a  point  outside  of  the  boundary  layer.  Using  the  requirement  for  constant  flow  rate 
to  eliminate  5,*  from  (3.20),  we  obtain,  after  a  little  rearranging, 

U"  ~  S*(l  +  etfj  +  e2K2  +  •••).  where  (32i) 

in,  -m,)  (K0  -1). 

Using  the  expansions  for  y*  and  0*  in  equation  (3.18),  we  obtain  on  collecting  powers  of  e  a 
series  of  problems  for  the  fn,  the  first  three  of  which  are 


2/"'  0, 

2/r  +fj"  -  /,7a'  -  v;'k  =  and 

2/r  +  /,/,'  -  2///;  +  3//7,  =  -2/2/2"  +  (/2')2  -  K\  -  2K2. 


(3.22) 
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Boundary  conditions  for  (3.22)  are  f„(0)  =  fn'(0)  =  0  and  fn'  -»  as  rj  — >  <»,  which  require  that 
t|/‘  satisfy  the  no  slip  condition  at  the  channel  wall  and  match  with  the  outer  flow  0*  at  the  upper 
edge  of  the  boundary  layer. 

The  first  of  equations  (3.22)  is  that  for  Blasius  flow  over  a  plate;  this  and  the  problems 
for  the  correction  terms  are  solved  in  a  straightforward  manner  using  a  multiple  shooting  routine 
(DBVPMS  from  IMSL).  The  base  flow  Uo  is  then  given  by 

u0(z-,x)~S(f{  (q)  +  e/2'(r| )  +e2//(Tl)  +  •”).  (3'23) 

where  by  definition,  e(x)  =  xm  /  (1^  S)  and  q  =  (z  +  1)  /  e(x).  In  the  derivation  of  the  Orr- 
Sommerfeld  system,  Uq  was  assumed  to  be  a  function  of  z  only;  we  thus  evaluate  e  for  a  fixed 
value  of  x  =  Xo  to  eliminate  the  dependence  on  the  axial  coordinate.  Because  e  varies  as  x01/2  it 
is  then  necessary  to  impose  a  downstream  limit  on  to  maintain  the  validity  of  the  perturbation 
solution  for  Uo.  Following  Schlichting  we  choose  to  restrict  Xo  so  that  e  is  sufficiently  small  that 
the  successive  terms  of  the  solution  for  Uq  remain  well  ordered,  that  is,  f/  >  ef2'  >  e^'.  This 
requires  that  e  be  less  than  0.1455,  or  by  extension  that  we  choose  Xo  <  (0.1455)2 S.  If 
additional  terms  are  included  in  the  expansion  for  Uo  the  downstream  constraint  on  Xq  becomes 
more  severe,  (e.g.,  Collins  and  Schowalter  1962  cite  e  <  0.0707  (Xo  <  0.005  R*  S)  for  a  seven 
term  expansion.) 

The  solution  (3.23),  normalized  by  the  magnitude  of  S,  is  shown  in  figure  3.2.  In  3.2a 
Uo  and  f,'  (the  Blasius  solution)  are  plotted  as  functions  of  z  and  compared  with  the  finite 
difference  solution  to  the  boundary  layer  equations,  for  a  representative  set  of  parameter  values. 
The  agreement  between  the  perturbation  and  numerical  solutions  is  seen  to  be  good.  In  3.2b  the 
values  of  f/,  e  f2'  and  e2  f3'  are  shown  for  the  same  choice  of  parameters. 

The  advantage  of  using  the  perturbation  solution  for  u^  (3.23)  arises  from  its  formulation 
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in  teims  of  the  Blasius  similarity  variable.  With  the  velocity  profile  tabulated  as  a  function  of 
this  independent  variable,  the  value  of  Uo  may  be  found  for  any  given  position  (x^z)  by 
calculating  the  value  of  rj  corresponding  to  that  point  (at  the  Reynolds  number  and  flow  speed 
being  considered)  and  then  interpolating  using  a  cubic  spline  routine  (from  Press  et.al.  1990)  to 
find  the  corresponding  value  for  Uq.  Solutions  to  the  stability  problem  (3.13),  (3.15-17)  with  the 
perturbation  solution  for  Uq  are  obtained  analytically  for  long  waves  and  numerically  for 
disturbances  with  arbitrary  wavelengths  in  Chapter  5. 


CHAPTER  4:  ASYMPTOTIC  AND  EXACT  SOLUTIONS  FOR 
A  PLUG  FLOW  BASE  STATE 

§4.1  Introduction  and  equations 

The  analytical  study  of  the  Orr-Sommerfeld  equation  for  non-uniform  base  velocity 
profiles  is  a  difficult  task  even  for  the  case  of  a  system  with  rigid  walls,  as  it  involves  the  solution 
of  a  non-constant  coefficient,  complex  valued,  fourth  order  eigenvalue  problem.  To  facilitate  the 
development  of  an  analytical  solution,  we  therefore  choose  to  approximate  the  developing  flow 
profile  with  a  uniform  (plug)  flow,  as  indicated  in  §3.3.  In  this  case,  it  is  possible  to  solve  the 
Orr-Sommerfeld  system  using  the  method  of  matched  asymptotic  expansions  to  obtain  an  explicit 
solution  for  the  eigenfunction  9  and  eigenvalue  c.  As  the  problem  with  plug  flow  has  constant 
coefficients,  it  is  also  possible  to  obtain  an  exact  solution,  which  determines  the  eigenfunction  and 
eigenvalues  implicitly.  We  present  the  matched  asymptotic  and  exact  solutions  in  §4.2  and  §4.3, 
respectively,  and  show  results  for  both  cases  in  §4.4. 

For  a  plug  base  flow  with  nondimensional  flow  speed  S,  the  Orr-Sommerfeld  equation 
(3.13)  becomes 

k(S  -  c)( 9"  -  k2 9)  +  _f_(9(,v)  -  2 k2<p"  +  *49)  =  0.  (4.1) 

Rw 

At  the  wall  (z  =  -1),  boundary  conditions  on  9  are  (3.15)  and  (3.16),  or 

(-(/fcc)2  -  2ikcG  +  BkA  +  Tk2  +  IH9  + 

(S  -  c)2 k(p'  +  (9'"  -  3 *V)  =  0  (4'2) 

Rw 

and 
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dk( p'  +  i£(< p"  +  *2<p)  =  o, 

Rw 

and  at  the  midline  (z  =  0)  <p  must  satisfy  symmetry  or  antisymmetry  conditions  (3.17), 


(4.3) 


<p  =  <p"  =  0  (symmetric),  or  <p'  =  <p"'  =  0  (antisymmetric) . 

In  the  equations  following,  we  denote  the  term  (-(k  c)2M  -  2ikcG  +  Bk4+  T  k2  +  1) 
appearing  in  (4.1)  as  A,  a  function  of  k  and  c. 

§4.2  Asymptotic  solution 

The  plug  base  flow  introduced  in  §4. 1  is  valid  in  the  limit  of  Reynolds  number  1^  — » 
so  that  it  is  appropriate  to  expand  (p  and  c  in  power  series  in  the  small  parameter  e  s  Rw1/2,  as 


9  -  <Po  +  e(Pi  +  e  92  + 

c  ~  c0  +  ec1  +  e2c2  +  •••. 

Plugging  (4.5)  into  (4.1)  and  letting  e  go  to  zero,  we  obtain  the  leading  order  problem 


(4.5) 


k(S  -  c0)(%"  -  k2%)  =  0  (4-6) 

with  midline  boundary  conditions  (4.4).  Because  in  the  limit  e  ->  0  the  two  highest  derivatives 
in  (4.1)  are  lost,  the  perturbation  is  singular  and  it  is  not  possible  for  this  reduced  problem  to 
satisfy  all  of  the  boundary  conditions.  A  boundary  layer  is  introduced  to  admit  the  effect  of 
viscosity  at  the  wall,  so  that  the  solution  to  the  outer  problem  must  satisfy  only  the  midline 
conditions;  the  inner  (boundary  layer)  solution  satisfies  the  conditions  at  the  wall  and  matches  to 
the  outer  solution  at  the  edge  of  the  boundary  layer.  Equation  (4.6)  is  the  Raleigh  equation  (see, 
e.g.,  Drazin  and  Reid  1981);  to  consider  inviscid  stability  we  would  require  that  <p  solve  this  and 
the  first  of  either  the  symmetric  or  antisymmetric  boundary  conditions  (4.4),  and  then  demand  that 
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the  normal  boundary  condition  at  the  wall  (4.2)  be  satisfied. 

The  solution  of  (4.6)  when  the  wave  number  k  is  non-zero  and  the  leading  order  phase 
speed  Co  is  not  equal  to  the  flow  speed  S  is,  up  to  an  arbitrary  multiplicative  constant, 

<p0  =  sinh(kz)  or  <p0  =  cosh(kz)  (4-7) 

(for  symmetric  or  antisymmetric  disturbances,  respectively).  If  the  wave  number  k  or  the 
difference  S  -  c  is  small  (0(e2)),  it  is  no  longer  possible  to  construct  an  asymptotic  solution 
differing  significantly  from  the  exact  solution  presented  in  §4.3,  and  we  thus  do  not  deal  with 
these  possibilities  here  ( c.f  §4.3).  The  solution  (4.7)  is  in  fact  the  exact  solution  to  (4.1)  for  the 
region  away  from  the  wall,  so  that  all  correction  terms  <pn  (n  >  0)  are  identically  zero. 

To  find  the  inner  solution,  a  new  independent  variable  £  is  introduced  to  stretch  the 
boundary  layer  near  the  wall.  With  the  small  parameter  5(e)  defining  the  boundary  layer  width, 
£  is  defined  by  £  =  (1  +  z)  /  5(e).  We  choose  5  to  retain  the  highest  derivative  term  in  the  Orr- 
Sommerfeld  equation  as  e  tends  to  zero,  which  requires  5(e)  =  e;  this  choice  in  fact  motivates 
the  use  of  e  as  our  expansion  parameter.  Rewriting  (4.1)  in  terms  of  the  inner  variable  £  and 
letting  0(0  =  cp(z),  we  obtain 

(<D<iv)  -  ik(S  -  c)0")  +  t2k\ik(S  -  c)0>  -  O")  +  e4*40  =  0. 

Similarly  rewriting  the  wall  boundary  conditions  (4.2)  and  (4.3)  gives 

(k(c  -  S)20'  -  i(c  -  5)0"')  +  eifcOA  +  3iV*2(c  -  S)<D'  =  0  <4-9) 

and 

dk&  +  eic<X>"  +  e3i/k2c< I>  =  0,  (4-10) 

both  at  £  =  0.  These  give  two  boundary  conditions  on  <J>;  the  inner  solution  must  also  match  with 
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the  outer  solution  <p  at  the  edge  of  the  boundary  layer.  This  matching  condition  is  formulated 
mathematically  below. 

Expanding  O  in  powers  of  e  as 


<I>  ~  <J>0  +  e<J>j  +  e202  +  •••  (4.11) 

and  introducing  this  and  the  expansion  for  c  (4.5)  into  the  inner  equation  and  boundary  conditions, 
a  series  of  problems  for  On  (n  =  0,1,...)  are  found  at  successive  orders  in  e.  At  leading  order, 

S£<D0  =  0,  where  $ N(S  -  c0)jll  -  ik(S  -  c0)2-l_,  (412) 

0  u  dz4  dz2 

with  boundary  conditions 

/?O0  =  0,  where  R  s  k(c0  -  S)JL  -  i^L—,  at  C  =  (4-13) 

1  dz  dz 

and 

dkO0'(0)  =  0 .  (4-14) 

Equation  (4.12)  is  homogeneous  with  constant  coefficients,  and  has  the  general  solution 

=  «01  +  a0zC  +  a03eP‘t  +  a04e^’  (4-15) 

where  p±  2  ±  ( i  k  (S  -  Co))w.  To  match  with  the  outer  solution  we  use  the  method  of 
intermediate  limits  (Kevorkian  and  Cole  1981)  and  introduce  the  intermediate  variable  %  = 
(1  +  z)  /  T|  (where  e  |ln  e  |  « tl(e)  «  1,  with  11(e)  ->  0  as  e  ->  0).  Rewriting  (4.15)  and  (4.7) 
in  terms  of  4  and  then  reexpanding  (4.7)  for  small  t|,  we  obtain 
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r,  p.Jli  P*J2i 

®,  ■  *  *•>'  '  *  “«e  ‘  “d  (4  ,6) 

%  -  r(*)  -  n*?Y'(*) *  *  •••. 

where  for  notational  convenience  the  term  y  (k)  =  -sinh(k)  or  cosh(k)  (for  symmetric  or 
antisymmetric  disturbances,  respectively)  has  been  introduced.  The  appropriate  matching 
condition  is  then 


Lint  (0-<p)=0,  (4.17) 

•-►o  v 

\  fixed 

which  requires  that  =  y  (k)  and  =  0.  The  tangential  boundary  condition  (4.14)  then 

demands  that  the  constant  a^  also  be  zero,  so  that 

<D0  =  y(*).  (4-!8) 

The  normal  boundary  condition  (4.13)  is  trivially  satisfied  by  (4.18),  and  the  eigenvalue  c0 
remains  undetermined;  the  leading  order  dispersion  relation  is  found  at  0(e). 

With  (4.18),  we  have  at  0(e) 

S£<X>,  =0,  (4-19) 

=  -A0*y,  (420) 

(where  Ag  =  A(k,c0)),  and 


<**0/(0)  =  0.  (4.21) 

Equation  (4. 19)  for  is  the  same  as  (4. 12)  for  O0,  so  that,  after  applying  the  first  order  matching 


condition 
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Lim  -1(0  -  <p)  =  0 

t~>0  0 

\  fixed 


and  the  tangential  boundary  condition  (4.21)  we  find 


(4.22) 


<t>,  ■  -ky'(k)^  +  e*-?.  (4-23) 

The  normal  boundary  condition  (4.20)  then  becomes  after  a  little  algebraic  manipulation  a 
quadratic  in  Cq  which  may  be  solved  to  give  the  leading  order  dispersion  relation 


S  -  iGP(k)  +  1 

\  +  MkP(k)  1  +  MkP(k) 


(S  -  iGP(k))2  - 


(1 


+  MkP(k)) (S2  -  (Bk*  +  Tk2  +  l)P(k) 
k 


(4.24) 


Here  P(k)  is  defined  by  P(k)  m  y(k)/y'(k)  (=  tanh(k)  or  cotanh(k)  for  symmetric  and 
antisymmetric  disturbances,  respectively).  As  the  limit  of  e  going  to  zero  corresponds  to  the  limit 
of  no  viscosity,  (4.24)  is  simply  the  inviscid  dispersion  relation  for  the  problem  (see,  e.g.  Grotberg 
and  Reiss  1984,  setting  their  viscous  correction  f  to  zero). 

To  find  the  first  (viscous)  correction  to  this  dispersion  relation  it  is  necessary  to  continue 
to  order  e2.  Using  the  expressions  for  O0  and  Oj,  (4.18)  and  (4.22),  and  collecting  terms  of 
0(e2),  <I>2  must  satisfy  the  inhomogeneous  equation 


se<D2  *  ik2cxy'^_e^  +  (2k2  -  ik2(S  -  c0))y 
with  boundary  conditions 

+  2k.2 Cj (M kc0  *  iG)y 


(4.25) 


(4.26) 


and 
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dk< D2'(0)  =  -ikc0  y'p_.  (4-27) 

A  particular  solution  for  (4.25)  is  found  using  the  method  of  undetermined  coefficients  to  be 

=  _  «*ac,Y'(*)  +  ±^£2  (kh  (4.28) 

2p  2p_2  2 

so  that,  after  matching  and  applying  (4.27)  as  at  earlier  orders, 


<*>,  =  ( 


ik,2cxy'{k)  ic0y'(k) 


)e**  -  +  h2^y(k).  (4.29) 


Inserting  (4.29)  into  the  normal  boundary  condition  (4.26)  gives  the  first  correction  to  the 
dispersion  relation. 


M(kc0)2  +  2ikc0G  -  ( Bk 4  +  T*2  +  1) 
2P_(S  -  (MkP(k)  +  l)c0  -  iGP(k)) 


Recalling  that  P_  =  -(/  k  (S  -  c0))1/2,  it  appears  at  first  sight  that  this  expression  is  singular  as  S 
approaches  c0.  However,  because  of  the  definition  of  c0  (4.24),  in  this  case  the  numerator  of 
(4.30)  goes  to  zero  more  rapidly  than  (S  -  Cg)1/2,  and  there  is  no  singularity.  A  singularity  may 
still  arise  if  the  remaining  term  in  the  denominator,  S  -  (M  k  P(k)  +  1)  Cq  - 1  G  P(k),  vanishes, 
which  for  G=0  (an  undamped  wall)  and  c0  strictly  real  occurs  as  the  flow  speed  S  approaches  the 
value  (M  k  P(k)  +  1)  c0.  The  leading  order  dispersion  relation  is  inviscid,  so  that  when  there  is 
no  wall  damping  Cq  will  in  fact  be  real  valued  until  the  onset  of  instability,  and  the  denominator 
may  thus  be  expected  to  go  to  zero  for  some  value  of  flow  speed.  When  wall  damping  is 
included,  as  it  will  be  in  all  physical  applications,  this  term  becomes  complex  valued,  so  that  it 
is  unlikely  that  it  will  completely  vanish.  For  small  wall  damping,  however,  the  imaginary  part 
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of  Co  is  small,  and  the  imaginary  part  of  this  term  will  also  be  small;  thus  for  lightly  damped  walls 
the  denominator  may  become  sufficiently  small  for  ct  to  become  0(e'‘)  and  the  asymptotic 
expansion  may  thus  break  down. 

Because  the  results  for  the  asymptotic  solution  developed  in  this  section  and  those  for  the 
exact  solution  in  §4.3  are  similar,  we  defer  further  discussion  of  the  dispersion  relation  (4.24, 
4.30)  until  §4.4. 

§4.3  Exact  solution 

For  plug  flow  the  Orr-Sommerfeld  equation  (4.1)  has  constant  coefficients,  so  that  an 
exact  solution  may  be  obtained  in  terms  of  exponentials.  Two  of  the  four  solutions  are  those  of 
the  inviscid  (outer)  solution  derived  above,  sinh(kz)  and  cosh(kz)  (which  are  appropriate  for 
symmetric  and  antisymmetric  disturbances,  respectively).  To  determine  the  remaining  two 
(viscous)  solutions,  let  <p  =  emkz;  inserting  this  into  (4.1)  gives 

k3(S  -  c)(m 2  -  1)  +  _L*4(m4  -  2 m2  +  1)  =  0.  (4.31) 

Rw 

Recognizing  that  m  =  ±1  correspond  to  the  inviscid  solutions  already  isolated,  and  assuming  k  * 
0  (we  return  to  this  possibility  below),  the  general  solution  for  <p  is  then  given  by 

i  R 

cp  ■  A^ikz)  +  A2g(kmz) ,  where  m2  -  1  +  __1(S  -  c),  (4.32) 

and  we  have  taken  g  (x)  =  sinh(x)  or  cosh(x)  according  to  whether  symmetric  or  antisymmetric 
disturbances  are  being  considered. 

For  simplicity  of  analysis,  we  assume  that  the  tangential  boundary  condition  at  the  wall 
may  be  replaced  by  that  corresponding  to  no  horizontal  motion,  (Uo(-l)  -  c)  cp'  -  Uo'(-l)  9  =  0  (c./. 
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§3.2),  as  this  assumption  makes  little  difference  to  the  solution  obtained.  (This  is  verified  with 
the  numerical  solution  of  the  full  system  formulated  in  chapter  5.)  Application  of  this  condition 
and  the  normal  wall  boundary  condition  (4.2)  yields  a  second  order  linear  homogeneous  system 
for  the  coefficients  Aj  and  A2, 


*(-*) A 
(5  -  c)g'(-k) 


g(-km)A  -  (5  -  c)2 kmg'(-km) 
(S  -  c)mg'(-km) 


0. 


(4.33) 


To  obtain  (4.33)  the  relation  g'"(kx)  =  k2g'(kx)  has  been  used.  The  eigenvalue  c  must  be 
chosen  to  zero  the  determinant  of  the  coefficient  matrix  on  the  left  hand  side,  which  gives  the 
dispersion  relation 


km(S  -  c)g'(-km)  g'(-k)  x 
/ 

A 


v,  t- 


P{-km)  _  P(-k)  _  (S  _  c)2 

km  k 


=  0. 


(4.34) 


In  (4.34),  as  in  the  asymptotic  solution,  P(x)  has  the  value  tanh(x)  or  cotanh(x)  for  symmetric  or 
antisymmetric  disturbances.  Solutions  of  this  dispersion  relation  arec  =  S,  c  =  S-ik  Re1,  and 
c  chosen  to  zero  the  parenthesized  term  in  the  relation.  The  first  two  of  these  possibilities  always 
yield  stability,  so  that  it  is  the  last  that  is  of  interest  Values  of  c  zeroing  the  parenthesized  term 
are  found  in  a  straightforward  manner  using  a  numerical  root  finder  (DNEQNF  from  IMSL). 
Note,  however,  that  c  =  S  is  also  a  solution  for  this  expression  (when  c  =  S,  m  =  1,  and  terms 


involving  P  cancel);  this  results  in  some  difficulty  in  the  numerical  resolution  of  other  roots  when 
they  approach  S. 

Information  about  the  limit  of  infinite  channel  width  for  (4.34)  may  be  obtained  by 
scaling  lengths  in  the  problem  on  L  =  Lw  =  (D*/E*  )1W,  instead  of  b*  ( c.f.  §3.1).  After 
appropriately  redefining  the  wall  and  fluid  parameters  (1^,  M,  etc:,  note  that  B=1  for  L=Lw),  the 
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dispersion  relation  (4.34)  is  changed  only  through  the  introduction  of  a  factor  of  (b*/Lw)  in  the 
arguments  of  g  and  P.  As  seen  in  §4.4,  symmetric  disturbances  are  the  least  stable,  and  we  thus 
consider  P(x)  =  tanh(x)  here.  In  the  limit  of  infinite  channel  width,  (b*/Lw)  -*  so  that 
P(kmb7Lw)  and  P(kb’/Lw)  both  go  to  1.  The  parenthesized  term  in  (4.34)  thus  becomes 

*  0,  (4-35) 

where  S  and  c  in  (4.35)  are  scaled  on  fl  =  (E'Lw/p  )1/2  and  k  on  Lw l.  Assuming  that  c  is  real 
(i.e.  that  points  being  considered  lie  on  the  neutral  stability  curve)  and  squaring  once  to  get  rid 
of  the  square  root  in  the  definition  of  m,  (4.35)  becomes  on  separating  real  and  imaginary  parts 


_L  +  1 

km  k 


-  (S  -  c)2 


k'q'  +  2RwA,q 2  -  2kARq  -  —JHARA,  =  0  and 
Rwkq*  -  2RwARq 2  -  2kA,q  +  Ol(A2R  -  A2,)  =  0. 


(4.36) 


In  (4.36)  q=S-candA  =  AR  +  i'A,.  Equations  (4.36)  are  solved  using  Newton’s  method  (we 
again  use  DNEQNF  from  IMSL);  results  are  shown  in  §4.4. 

If  the  wave  number  k  goes  to  zero,  the  Orr-Sommerfeld  equation  degenerates  to  <p°v)  = 
0,  and  the  boundary  conditions  at  the  wall  to  i  Re'1  (S  -  c)  <p"'  =  0  and  (S  -  c)  <p'  =  0.  For 
symmetric  disturbances  the  only  non-trivial  solution  is  c  =  S,  which  is  stable.  For  antisymmetric 
disturbances  the  only  eigenfunction  is  trivial,  and  hence  stable;  this  corresponds  to  the  movement 
of  the  fluid  across  the  midline  of  the  channel  as  a  solid  body.  Thus,  unlike  the  case  of  a 
developing  flow  profile  (c.f.  Chapter  5)  the  limit  k  ->  0  is  not  of  interest  for  the  purposes  of 
determining  the  onset  of  instability. 


40 


§4.4  Results 

The  asymptotic  dispersion  relation  (4.24),  (4.30)  has  two  roots,  the  flutter  modes  of  the 
system.  At  zero  flow  speed  these  have  equal  growth  rates  and  phase  speeds  of  equal  magnitude 
but  opposite  sign;  the  root  with  positive  phase  speed  we  call  root  1  and  the  other  root  2.  As  the 
flow  speed  is  increased,  root  2  is  slowed  until  a  flow  speed  SD,  at  which  it  reverses  direction  and 
travels  downstream.  In  an  undamped  system  (when  wall  damping  G  and  viscosity  e  are  0)  the 
phase  speeds  of  roots  1  and  2  coalesce  at  a  higher  flow  speed,  SF,  at  which  point  the  growth  rates 
of  the  roots  split  from  zero  to  become  positive  and  negative,  heralding  instability;  this  transition 
is  shown  in  figure  4.1.  As  noted  in  Chapter  2,  the  addition  of  wall  damping  breaks  the  unstable 
equilibrium  existing  for  flow  speeds  S  between  SD  and  SF,  so  that  instability  appears  for  root  2, 
at  S  =  SD.  The  addition  of  viscosity  may  restabilize  the  system  so  that  flutter  reappears,  at  a  flow 
speed  larger  than  SD,  as  shown  by  Grotberg  and  Reiss  (1982,  1984);  this  is  shown  in  figure  4.2, 
which  shows  root  2  of  the  asymptotic  and  exact  dispersion  relations. 

The  agreement  between  the  asymptotic  and  exact  solutions  is  clearly  very  good  in  figure 
4.2,  but  this  changes  when  the  wall  damping  is  decreased  slightly,  as  in  figure  4.3  (which  shows 
roots  1  and  2  of  the  dispersion  relation  in  4.3A  and  B,  respectively).  For  this  case  the  exact 
solution  shows  that  the  wall  damping  is  no  longer  sufficient  to  destabilize  root  2  (which  for  S 
between  SD  and  SF  is  Type  A,  as  defined  by  Benjamin  1963  and  Landahl  1962,  and  discussed  in 
Chapter  2),  so  that  root  1  (which  is  Type  B)  is  instead  that  which  goes  unstable.  The  asymptotic 
solution,  however,  still  predicts  stability  for  root  1  and  instability  for  root  2.  Figure  4.4  shows 
the  more  physical  case  of  lower  wall  damping,  for  which  the  singularity  anticipated  in  the 
asymptotic  solution  is  seen.  This  singularity  may  be  expected  any  time  the  wall  damping  is  small; 
however,  in  the  cases  in  which  we  are  interested,  air  flow  in  flexible  tubes  and  the  lung  airways, 
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Figure  4.2:  Plug  flow  dispersion  relation  for  damped  system,  root  2,  showing 
exact  (solid  curves)  and  asymptotic  (dashed  curves)  solutions.  •  real(c); 
o  imag(c).  M=3190,  B=542,  T=0,  G=16,  Rw=2230,  k=0.1 


O  O 


0.5 


1.0 

Flow  Speed 

Figure  4.3:  Plug  flow  dispersion  relation,  showing  difference  between  exact 

(solid  curves)  and  asymptotic  (dashed  curves)  solutions.  •  real(c); 
imag(c).  4. 3 A  shows  root  1,  4.3B  root  2.  M=3190,  B=540,  T=0, 
=15, 1^=2230,  k=0.1. 


•  real(c), 
o  imag(c) 


Figure  4.4:  Plug  flow  dispersion  relation,  showing  singularity  in  asymptotic 

solution.  Solid  curve  gives  exact  solution,  dashed,  asymptotic. 

•  real(c);  o  imag(c).  Root  1.  M=3190,  B=542,  T=0,  G=0.054, 

Rw=2230,  k=0.1. 
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the  solution  remains  valid  for  flow  speeds  up  to  and  including  the  onset  of  instability. 

The  discrepancy  between  the  behavior  of  the  exact  and  asymptotic  solutions  of  the 
problem  as  wall  damping  is  decreased  is  a  direct  result  of  the  assumptions  used  in  the  formulation 
of  the  latter.  As  the  expansion  parameter  e  in  the  asymptotic  solution  goes  to  zero  the  problem 
becomes  inviscid,  so  that  at  leading  order  the  solution  is  that  for  inviscid  flow.  Implicit  in  this 
expansion,  therefore,  is  the  assumption  that  the  effects  of  viscosity  may  be  included  as  a  small 
correction  to  the  inviscid  result,  an  assumption  that  is  clearly  not  true  in  figure  4.3.  Thus  as  long 
as  wall  damping  is  sufficiently  large  (by  comparison  with  viscosity)  to  maintain  the  stability  of 
root  1  (the  Type  B  flutter  mode)  the  asymptotics  agree  well  with  the  exact  solution.  However, 
when  the  wall  damping  is  smaller,  so  that  root  1  is  unstable,  there  is  an  order  one  discrepancy 
between  the  viscous  and  inviscid  solutions  for  larger  flow  speeds,  and  the  asymptotics  may  no 
longer  be  able  to  accurately  model  the  system.  Nonetheless,  for  walls  with  sufficiently  small 
damping  and  large  mass  ratios  (as  in  figure  4.4  and  noted  above),  the  asymptotic  solution  is  able 
to  accurately  predict  instability  for  root  1  before  the  appearance  of  the  singularity,  and  before  the 
predicted  behavior  for  roots  1  and  2  differs  significantly  from  that  in  the  exact  solution.  For  cases 
of  interest,  root  1  becomes  unstable  as  shown  in  figure  4.4,  before  instability  is  seen  for  root  2. 

In  figure  4.5  we  show  curves  of  neutral  stability  (lines  along  which  the  growth  rate  of  the 
disturbances  vanish)  for  the  plug  flow  base  state,  for  symmetric  (4.5A)  and  antisymmetric 
instabilities  (4.5B).  These  give  the  flow  speed  at  which  instability  appears  as  a  function  of  wave 
number,  points  below  (above)  the  curve  are  stable  (unstable).  The  minima  of  the  curves  give  the 
flow  speed  SCRF  at  which  flutter  appears,  and  the  corresponding  critical  wave  number  kCR. 
Comparing  4.5A  and  B  shows  that  the  symmetric  mode  is  the  least  stable,  which  is  generally  true. 
In  both  4.5A  and  B  curves  for  both  the  exact  and  asymptotic  solutions  are  plotted,  and  give 
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similar  but  not  identical  results;  the  agreement  between  the  two  is,  however,  better  for  the 
symmetric  mode  than  the  antisymmetric  mode.  In  that  the  symmetric  disturbances  become 
unstable  before  antisymmetric  ones  the  larger  discrepancy  in  the  latter  case  is  physically 
immaterial.  Both  symmetric  and  antisymmetric  disturbances  are  seen  to  be  stable  for  long  wave 
lengths  (small  k),  as  noted  in  §4.3. 

Several  neutral  stability  curves,  for  different  values  of  the  wall  damping  G  and  mass  ratio 
M,  are  shown  in  figure  4.6.  Note  that  because  the  mass  ratio  appears  in  the  definition  of  G,  to 
isolate  the  effect  of  variation  of  M  it  is  in  this  case  necessary  to  allow  G  to  change  while  holding 
6  (=  b*  G*  /  (1)  fixed.  As  the  symmetric  mode  is  unstable  first,  these  curves  are  for  the  symmetric 
instability.  For  both  values  of  the  mass  ratio  shown  (indicated  by  filled  and  open  symbols, 
respectively),  increasing  the  wall  damping  is  stabilizing,  which  may  be  expected  as  the  instability 
is  for  these  parameter  values  a  Type  B  instability,  which  is  stabilized  by  wall  damping. 
Increasing  6  also  increases  the  wavelength  of  (decreases  k  for)  the  instability.  For  the  larger  mass 
ratio  (filled  symbols),  the  instability  moves  to  still  smaller  k.  For  smaller  values  of  the  wall 
damping,  increasing  the  mass  ratio  is  destabilizing,  which  is  consistent  with  the  fact  that  the  flutter 
instability  is  destabilized  by  wall  inertia.  As  both  the  wall  damping  and  M  are  increased, 
however,  the  wavelength  of  the  instability  becomes  long  enough  that  the  system  is  stabilized  by 
the  band  of  long-wave  stability  noted  above  and  in  §4.3;  thus  for  sufficiently  large  wall  damping 
increasing  the  wall  inertia  is  in  fact  stabilizing.  This  is  shown  in  figure  4.7  by  plotting  just  the 
minima  of  the  neutral  stability  curves,  SCRF ,  for  different  values  of  the  wall  damping  G.  For  the 
smaller  (more  physically  relevant)  values  of  the  wall  damping  increasing  M  is  destabilizing,  but 
this  is  not  the  case  for  more  heavily  damped  walls.  This  figure  also  compares  the  results  for  the 
exact  and  asymptotic  solutions;  the  two  are  seen  to  agree  reasonably  well,  with  the  same  essential 
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Figure  4.6:  Plug  flow  neutral  stability  curves  for  different  M,  G,  showing  long 
wave  stabilization.  Open  symbols  give  M=500,  filled,  M=5000; 
o  6x5000=0.1567;  v  6x5000=1.567;  □  6x5000=4.702.  Symbols  are 
located  at  the  minima  of  the  neutral  stability  curves.  B=542,  T=0, 
1^=2230.  (6  s  b*  G*  /  u) 


Wall  Damping  5000  G 


Figure  4.7:  Plug  flow  critical  flow  speed  SCRF  (4.7 A)  and  flutter  frequency  (4.7B) 
as  functions  of  wall  damping,  o  M=500,  exact  solution;  •  M=5000, 
exact  solution.  Dashed  and  dotted  curves  give  asymptotic  solution  for 
M=500  and  M=5000,  respectively.  B=542,  T=0,  1^=2230. 
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features  being  present  in  either  case,  though  the  values  of  SCRF  differ  slightly.  As  a  result, 
subsequent  figures  show  only  the  exact  solution,  with  the  understanding  that  the  asymptotic  result 
may  be  expected  to  retain  the  same  features  for  slightly  different  critical  flow  speeds.  Figure  4.7B 
shows  how  the  variation  of  G  and  M  affects  the  flutter  frequency  at  the  onset  of  instability; 
increasing  either  of  the  wall  mass  or  the  wall  damping  decreases  the  frequency,  which  in  either 
case  is  consistent  with  physical  intuition  and  observation.  This  also  shows  that  the  predictions 
for  the  critical  flutter  frequency  made  by  the  asymptotic  and  exact  solutions  are  in  very  close 
agreement.  Thus  if  the  flutter  frequency  is  the  object  of  interest,  the  asymptotic  solution  (which 
is  explicit  and  hence  easier  to  use)  may  be  used  with  little  or  no  loss  of  accuracy. 

Figure  4.8  shows  how  variation  of  the  wall  elastance  and  bending  stiffness  B  changes  the 
stability  of  the  system,  by  plotting  (in  4. 8 A)  SCRF  as  a  function  of  the  wall  elastance  E 
(nondimensionalized  on  (p*  v*2  /  b'^xlO4)  for  different  values  of  B.  For  this  figure  SCRF  is 
scaled  on  (v*  /  b')xl02,  and  the  frequency  (shown  in  4.8B)  on  v*/b*2.  Increasing  E  or  B 
stabilizes  the  flutter  instability,  and  increases  the  frequency.  The  frequency  increases 
approximately  as  E1/2,  which  is  consistent  with  the  results  of  Grotberg  and  Reiss  (1984).  In  that 
the  flutter  instability  depends  fundamentally  on  wall  elastance,  the  increase  in  SCRF  with  increasing 
E  is  as  expected;  in  the  limit  of  infinite  E  (a  rigid  wall),  the  flutter  instability  ceases  to  exist. 

In  figure  4.9,  we  show  the  effects  of  variation  in  the  half  channel  width  b,  for  different 
values  of  the  wall  elastance  E.  For  later  comparison  with  the  results  presented  in  Chapter  5,  we 
scale  velocity  on  (v*  /  x^xlO2,  length  on  x<,*  and  frequency  on  (v*  /  x*,*2),  where  Xq*  is  a  fixed 
distance  from  the  mouth  of  the  channel.  (We  have  chosen  x0*  to  be  equal  to  the  width  of  the 
channels  considered  in  figures  4.1-8.)  As  the  half  channel  width  is  increased,  the  effective 
viscosity  in  the  channel  decreases,  which  stabilizes  the  flutter  instability,  as  seen  in  figure  4.9 A. 


Wall  Elastance  E 

Figure  4.8:  Plug  flow  SCRF  (4.8A)  and  flutter  frequency  (4.8B)  as  functions  of 
wall  elastance  for  B=2.688xl04  (o)  and  B=2.688xl05  (•).  E*  scaled  on 
(p*  v*  /b‘3)xlO\  flow  speed  on  (v*/b*)xl02,  and  frequency  on 
(v*  /  b* 2 ).  Other  parameters  as  figure  4.4,  modified  for  scaling. 
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Figure  4.9:  Plug  flow  SCRF  (4.9A)  and  flutter  frequency  (4.9B)  as  functions  of 
half  channel  width  b*  /  Xq*  for  E*  /  (p*  v* 2  /  Xq*  3  )xlO*  =  500  ( v)  and  E  = 
900  (▼).  Velocity  and  frequency  scaled  on  (v*  /  Xq*  )xl02  and  (v*  /  Xq*2). 
Xq*  chosen  =  b*  in  figure  4.4;  other  parameters  as  figure  4.4,  modified  for 
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In  the  limit  of  very  wide  channels  the  critical  flow  speed  levels  off,  as  seen  at  the  end  of  the 
channel  widths  shown.  Increasing  the  wall  elastance  is,  as  seen  in  figure  4.8,  stabilizing.  In 
figure  4.9B  the  frequency  of  the  flutter  instability  at  onset  is  shown  as  a  function  of  the  half 
channel  width.  This  shows  that  as  the  half  channel  width  is  increased  through  moderate  values, 
the  flutter  frequency  decreases,  but  that  for  very  large  channel  widths  it  increases  again.  This 
change  in  behavior  may  be  understood  by  considering  the  neutral  stability  curves  giving  figure 
4.9,  which  are  shown  in  figure  4.10.  These  show  curves  for  increasing  b,  along  with  the  infinite 
channel  limit  (obtained  from  equations  (4.36)).  As  b  is  increased,  the  neutral  stability  curves 
agree  with  the  infinite  b  limit  for  disturbances  with  longer  wavelengths  (smaller  wave  numbers). 
Thus  the  critical  wave  number  kcR,  which  initially  decreases  with  increasing  b,  reverses  to  increase 
towards  the  infinite  channel  limit  for  very  wide  channels.  This  increase  results  in  the  increase  in 
flutter  frequency  shown  in  figure  4.9B.  There  is,  of  course,  always  a  difference  between  the  finite 
and  infinite  channel  width  curves  for  very  long  waves.  Note  that  the  channel  widths  relevant  to 
physical  tubes  (e.g.  those  of  the  experiments  of  Gavriely  et.al.  1989)  are  b  <  4  in  figure  4.9. 

As  we  are  specifically  interested  in  pulmonary  applications,  it  is  useful  at  this  juncture 
to  consider  the  effects  we  have  discussed  above  specifically  in  this  context.  Wheezing  lung 
sounds  are  thought  to  be  symptomatic  of  airway  flutter  (Grotberg  and  Davis  1980);  thus  the 
characteristics  that  distinguish  the  lungs  of  wheezing  patients  from  healthy  subjects  should  relate 
to  those  effects  that  destabilize  the  flutter  instability.  We  saw  above  that  these  destabilizing 
effects  include  reduction  of  the  wall  elastance  and  bending  stiffness,  and  reduction  of  the  half 
channel  width.  All  of  these  do  correspond  to  the  characteristics  of  asthmatic  or  lung  obstructed 
patients’  airways,  which  have  more  flexible  walls  and  are  narrower  than  those  of  healthy  subjects. 
We  directly  compare  our  theoretical  results  with  experimental  observations  in  Chapter  6. 
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Figure  4.10:  Plug  flow  neutral  stability  curves  for  different  half  channel  widths. 
Dotted  through  solid  curves  are,  respectively  b  =  2,  6.4,  8.6,  11.4  and 
infinite  width.  Scales  and  parameters  as  figure  4.9  (E*  /  (p*  v* 2  /  x0* 3 )  = 
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CHAPTER  5:  NUMERICAL  AND  ANALYTICAL  SOLUTIONS  FOR 
DEVELOPING  BASE  FLOW 

§5.1  Introduction  and  equations 

In  this  chapter  we  develop  solutions  to  the  Orr-Sommerfeld  system  derived  in  Chapter 

3  with  a  developing  base  flow.  An  analytical  solution  is  obtained  for  long  wave  disturbances  in 
§5.2,  and  a  numerical  solution  for  disturbances  of  arbitrary  wavelengths  is  formulated  and  tested 
in  §5.3  and  §5.4.  Results  are  presented  and  compared  with  the  plug  flow  solutions  of  Chapter 

4  in  §5.5,  and  a  short  discussion  of  the  limit  of  infinite  channel  width,  which  turns  out  to  be  of 
mathematical  interest,  appears  in  §5.6.  In  Chapter  6  the  solutions  presented  in  this  chapter  are 
compared  with  experimental  results.  Both  the  analytical  and  numerical  solutions  presented  below 
are  general  to  flow  profiles  other  than  the  developing  flow,  and  we  are  therefore  able  to  check 
our  results  through  consideration  of  problems  for  which  results  are  already  known. 

For  flow  profiles  Uq(z)  satisfying  the  no  slip  boundary  conditions,  such  as  the  developing 
flow  solution  (3.23),  the  Orr-Sommerfeld  equation  is 

k(u0  -  c)(cp"  -  *2(p)  -  *u0"<p  +  J_(<p'v  -  2*2(p"  +  *4cp)  =  0  (5.1) 

Rw 

and  the  boundary  conditions  (3.15)  and  (3.16)  are 

( ~(kc)2M  -  2i*cG  +  Bk*  +  Tk2  +  1)*9  + 

*c(c<p'  +  <<p)  -  il(< p'"  -  3/k2<p')  =  0  (52) 

Rw 

and 
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dk(cq>'  +  w0'q>)  +  iiL(<p"  +  *2<p)  =  0,  (5.3) 

R 


both  at  z  =  -1.  The  midline  boundary  conditions  (3.17)  remain 


<p(  0 )  =  <p"(0)  =0  (symmetric),  or  (5.4) 

<p'(0)  =  <p'"(0)  =  0  (antisymmetric). 


§5.2  Analytical  solution  for  long  waves 

While  it  is  in  general  difficult  to  solve  the  Orr-Sommerfeld  system  analytically,  as 
indicated  in  Chapter  4,  a  solution  for  the  case  of  long  waves  (k  -»  0)  is  more  tractable,  and  is 
developed  in  this  section.  We  consider  the  limit  of  no  horizontal  wall  motion  here  (i.e.,  reduce 
(5.3)  to  c<p'  +  Uo'q>  =  0  at  z  =  -1)  to  simplify  this  analysis;  comparison  with  the  numerical  solution 
(described  in  §5.3)  shows  that  the  results  obtained  are  quantitatively  the  same  as  those  obtained 
with  the  full  boundary  condition.  In  the  following  we  do  not  a  priori  specify  the  base  flow 
profile  tv,  specific  cases  are  considered  after  the  derivation.  For  k  «  1 ,  let 

<p  -  <p0  +  k  +  •  •  •  and  c  ~  c0  +  k  +  •  •  •  .  (5-5) 

Using  these  in  the  Orr-Sommerfeld  equation  (5.1)  and  letting  k  — »  0,  we  obtain 

<Po  ®  0,  (5-6) 

so  that  %  is  a  cubic  polynomial  in  z.  Imposing  the  boundary  conditions  (5.4)  and  the  leading 
order  forms  of  (5.2)  and  (5.3),  we  find 

<p0  ■  z  and  c0  *  u0'(z  =  -1 ) .  (5-7) 

Here  we  have  normalized  %  so  that  |<p0(-l)|  =  1.  The  solution  (5.7)  is  for  symmetric 
disturbances;  we  address  the  antisymmetric  mode  below.  We  see  from  (5.7)  that  c0  is  strictly  real. 
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and  so  must  continue  to  the  next  order  in  k  to  evaluate  the  stability  of  the  system. 
At  order  k,  the  Orr-Sommerfeld  equation  (5.1)  is 


(5.8) 


9i  -  -<***« o 

and  (p!  is  thus  the  sum  of  another  cubic  polynomial  in  z  and  a  particular  solution  satisfying 

%r  - 1*.(%  -  *<>  <5-9) 

(after  integration  by  parts  once;  further  integration  is  dependent  on  the  choice  of  the  base  flow 
u0).  Applying  the  boundary  conditions  (5.4)  and  the  order  k  terms  of  (5.2),  we  find 


9,  =  -9„>(z=0)  +  Bi 


so  that,  from  (5.3), 


2  6  <(z  — 1) 


iR 


«0'(z=-l)  [  O^Cz-O)  -  |®,/(z-0)  -  ^,(2  —  1)  - 


«,;(*—!)  i  - 


3<(z  — 1) 


[1  -  (u0'(z  =  -l))2] 


(5.11) 

A 


where  we  have  written  (p1P  =  i  R*  <b1P  to  make  the  imaginary  character  of  C!  explicit. 

Equation  (5.11)  for  C!  is  strictly  imaginary,  so  that  cx  is  the  growth  rate  for  the 
disturbance.  To  find  the  critical  flow  speed  at  which  c,  becomes  positive,  it  is  necessary  to 
integrate  the  expression  (5.9)  for  tp1P.  We  consider  two  cases.  For  Poiseuille  flow,  =  1.5  S  (1  - 
z2),  and  we  may  integrate  to  find  <p1P  and  evaluate  (5.11)  analytically,  obtaining  c:  = 
0.067  i  R*  (18  S2  -  5).  Instability  appears  at  SCRLW  =  0.527;  at  this  point  Cq  =  Uo'(-l)  =  1.581. 
These  are  in  good  agreement  with  the  numerical  solution  (described  below),  which  gives  SCRLW  = 
0.525  and  c  =  1.575  (at  k  =  0.001).  For  developing  flow  (3.23)  we  numerically  integrate  (5.9) 
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using  Simpson’s  rule,  and  find  (for  1^=2230  and  x=l)  SCRLW  =  0.1387  and  Co  =  1.094;  again, 
these  agree  well  with  the  values  from  the  numerical  solution  of  the  full  problem,  SCRLW  =  0.1388 
and  c  =  1.093  (at  k  =  0.001).  The  existence  of  the  long  wave  instability  is  also  useful  for  the  full 
numerical  solution,  as  it  gives  starting  points  from  which  the  stability  calculation  may  proceed. 

The  effect  of  variation  of  the  half  channel  width  and  wall  elastance  on  the  long  wave 
stability  boundary  for  developing  flow  is  shown  in  figures  5.1  and  5.2,  which  show  the  long  wave 
critical  flow  speed  SCRLW  and  volumetric  flow  rate,  b  SCRLW,  as  functions  of  the  half  channel  width 
and  wall  elastance.  In  both  figures  flow  speeds  and  half  channel  width  are  scaled  on  (v  /  x*,* ) 
and  x o*,  respectively,  to  isolate  the  effect  of  changing  b;  by  showing  curves  for  increasing 
Reynolds  numbers  R*  (defined  as  t  Xq*  /  v*,  where  fi  =  (E*  Xo*  /  p')w),  the  effect  of  increasing 
elastance  is  also  seen.  These  figures  show  that  the  critical  flow  speed  increases  as  channel  width 
is  decreased,  but  that  the  critical  volumetric  flow  rate  decreases.  Increasing  the  wall  elastance 
increases  both  the  critical  flow  speed  and  flow  rate,  thus  stabilizing  the  system  significantly. 

Note  that  the  stability  bounds  derived  above  do  not  depend  on  the  wall  properties  other 
than  the  elastance;  the  effects  of  the  mass  ratio  M,  bending  stiffness  B,  etc.,  appear  in  the  stability 
calculation  only  at  higher  orders  in  k.  As  a  result,  variation  of  these  wall  parameters  may 
stabilize  the  system  only  to  the  point  at  which  the  long  wave  instability  is  critical.  Furthermore, 
it  is  possible  to  find  wall  parameters  such  that  the  long  wave  and  finite  wavelength  instabilities 
become  unstable  at  the  same  flow  speed,  resulting  in  a  codimension  two  bifurcation  point.  This 
is  shown  in  figure  5.1  by  also  showing  the  critical  flow  speed  for  flutter  (the  finite  wavelength 
instability),  SCRF,  as  a  function  of  b;  the  codimension  two  point  occurs  when  this  curve  intersects 
that  for  the  long  wave  instability.  Variation  of  other  system  parameters  changes  SCRLW  in  the 
same  manner  as  SCRF  (c./.  §5.5). 
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Figure  5.1:  Critical  flow  speeds  as  functions  of  half  channel  width  and  elastance. 
b*  scaled  on  x<,*,  S*  on  (v*  /  V  )xl02.  v,  ▼:  SCRLW;  o,  •: 
Scrf.  Open  symbols  give  =  2230,  filled  R^,  =  3000  (where  R*,  = 
(il  Xq*  /  v)  and  fl  =  (E*  x„*  /  p*)1/2).  Figure  5. IB  shows  the  frequency 
(scaled  on  (v'/xq*2)  of  the  instability  at  onset.  Other  parameters  as 
figure  5.7  (M  =  3190),  with  modified  scaling. 
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Figure  5.2:  Critical  volumetric  flow  rate  for  long  wave  instability  as  a  function 
of  half  channel  width  and  elastance.  Scales  and  parameter  values  as  in 
figure  5.1.  v  R*  =  2230;  *  Rw  =  3000;  □  R*  =  5000. 
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The  above  treatment  was  strictly  for  the  case  of  symmetric  disturbances.  When 
antisymmetric  disturbances  are  considered  the  analysis  proceeds  in  much  the  same  manner,  with 
the  midline  boundary  conditions  replaced  by  the  antisymmetric  forms  of  (5.4).  In  this  case  the 
leading  order  term  of  the  wall  condition  (5.2)  and  the  second  midline  condition  are  redundant, 
so  that  the  solution  for  the  leading  order  antisymmetric  eigenfunction  and  eigenvalue  <pA0  and  cA0 
becomes  a  one  parameter  family  of  solutions, 


%o  =  aao  +  daoz1  •  and  c*o  =  (^°  „  —  uo  (z  =  "I )  • 


2D 


(5.12) 


AO 


(for  normalization  we  may  specify  one  of  Aa0  or  DA0;  when  possible  we  take  Aa0  =  1)  provided 
the  free  parameter  DA0  is  non- zero.  If  DA0  is  zero  the  second  wall  condition,  (5.3),  requires  Aa0  = 
0,  and  the  solution  for  cpA0  is  trivial.  At  0(k),  condition  (5.2)  becomes  a  condition  on  DA0,  and 
requires  it  to  be  zero.  Thus  the  only  eigenfunction  is  the  trivial  one,  and  there  is  no  long  wave 
antisymmetric  instability.  This  is  similar  to  the  antisymmetric  long  wave  limit  for  plug  flow 
(§4.3);  in  the  long  wave  limit  the  fluid  may  be  moved  as  a  solid  body  across  the  midline, 
resulting  in  a  zero  disturbance  flow  field. 

It  may  be  noted  at  this  juncture  that  in  the  limit  of  long  waves  the  assumption  of  locally 
parallel  flow  made  in  the  derivation  of  the  Orr-Sommerfeld  equation  is  violated.  Thus  in  this 
limit  the  developing  flow  result  would  be  somewhat  different  if  non-parallel  effects  were 
included;  possible  differences  are  considered  in  the  discussion  of  Chapter  7.  Clearly  this  is 
irrelevant  for  Poiseuille  flow,  and  the  similarity  between  the  Poiseuille  and  developing  flow 
stability  results  suggests  that  the  developing  flow  result  is  at  least  qualitatively  accurate.  Note 
that  if  the  wall  boundary  conditions  (5.2,  5.3)  are  replaced  by  those  for  a  rigid  wall,  <p(z=-l)  = 
cp'(z=-l)  =  0,  the  solution  for  either  symmetric  or  antisymmetric  distuibances  is  trivial,  confirming 
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that  the  long  wave  limit  does  not  correspond  to  the  Tollmien-Schlichting  instability. 

§5.3  Numerical  solution  method 

Because  the  Orr-Sommerfeld  equation  is  notoriously  stiff,  standard  linear  shooting 
techniques  that  might  otherwise  be  used  to  obtain  a  solution  to  the  eigenvalue  problem  may  fail, 
especially  for  high  Reynolds  numbers  flows.  We  therefore  develop  in  this  section  a  multiple 
shooting  algorithm  with  orthonormalization  to  solve  (5.1)  with  boundary  conditions  (5.2-4) 
numerically.  This  orthonormalization  method  is  similar  to  that  proposed  by  Davey  (1973),  with 
modifications  to  allow  for  the  more  complicated  boundary  conditions  introduced  by  the  compliant 
wall. 

Linear  shooting  methods  solve  a  boundary  value  problem  as  an  initial  value  problem, 
using  supeiposition  of  solutions  to  satisfy  the  boundary  conditions  at  the  other  end  of  the  domain. 
To  accomplish  this  numerically,  (5.1)  is  written  as  a  first  order  matrix  equation  for  the  vector  p  = 
(<p  <p'  <p"  q>"')T,  as  p'  =  D  p,  where  the  entries  of  D  are  found  in  the  usual  manner  from  (5.1). 
Boundary  conditions  (5.2)  and  (5.3)  are  written  as  vector  conditions  on  p,  B,-  p  =  0  and  B2-  p  = 
0,  where  the  k*  components  of  B,  and  B2  are  the  coefficients  of  the  (k-l)fl  derivatives  of  q>  in 
(5.2)  and  (5.3),  respectively.  An  ’initial  condition’,  pM  for  p  is  chosen  to  satisfy  the  midline 
boundary  conditions  (5.4)  and  integrated  numerically,  according  to  the  vector  form  of  (5.1),  across 
the  channel  to  obtain  p  at  the  wall,  pw.  (We  use  the  variable  step  Runge-Kutta  routine  DIVPRK 
from  IMSL  for  the  numerical  integration.)  Because  the  problem  is  linear,  this  integration  may 
be  written  as  a  matrix  operation,  as  pw  =  R  p^  where  R  is  the  ’transfer  matrix’  for  the  problem. 
Note  that  the  jth  column  of  R,  rjt  is  found  by  integrating  the  jth  basis  vector  (where  = 
(1  0  0  0)T,  etc.)  across  the  interval.  The  boundary  conditions  (5.2)  and  (5.3),  in  vector  form,  then 
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Substituting  r,  for  the  jth  columns  of  R  and  plugging  in  the  symmetric  initial  condition  p„  = 
Aj  %  +  A2  e4,  (5.13)  becomes 
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For  a  non-trivial  solution,  then,  the  determinant  condition 


(5.14) 


(Brr2)(B^r4)  -  (B1T4)(B2-r2)  =  0  (5-15) 

must  be  satisfied.  This  condition  determines  the  eigenvalue  c,  which  is  iterated  from  an  initial 
guess  by  a  root  finding  algorithm  until  (5.15)  holds.  Clearly  the  same  procedure  may  be  applied 
for  the  antisymmetric  case,  replacing  the  basis  vectors  used  in  the  initial  condition  with  e,  and 
ej- 

This  procedure  will  woric  unless  the  vectors  r2  and  r4  in  (5.15)  are  parallel  or  nearly 
parallel,  that  is,  unless  the  integration  of  e2  and  e4  across  the  interval  yields  nearly  the  same 
solution.  To  see  how  this  may  occur,  consider  the  simple  system  <p"  -  m2  <p  =  0  with  m  large. 
The  general  solutions  for  9  are  e1”1,  and  any  initial  condition  will  be  a  linear  combination  of 
these.  Thus  as  any  general  initial  condition  (that  does  not  coincide  exactly  with  the  e  mz  solution) 
is  integrated  through  positive  values  of  z  it  will  rapidly  become  dominated  by  the  growing 
exponential,  and  solutions  for  differing  initial  conditions  will  thus  soon  become  almost  equal. 
The  Orr-Sommerfeld  system  similarly  has  exponentially  growing  and  decaying  solutions,  so  that 
even  though  the  initial  conditions  for  the  columns  in  the  transfer  matrix  R  are  orthonormal,  the 
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columns  themselves  may  still  be  almost  parallel.  Furthermore,  even  if  the  integration  of  the 
different  columns  of  the  transfer  matrix  proceeds  with  high  accuracy,  the  exponential  growth  will 
result  in  the  magnitude  of  the  columns  becoming  very  large,  so  that  the  numerical  roundoff  error 
in  the  evaluation  of  the  determinant  may  be  significant.  Systems  exhibiting  this  type  of  behavior 
are  described  as  stiff,  and  it  is  to  address  stiffness  that  orthonormalization  is  introduced. 

Orthonormalization  modifies  the  transfer  matrix  to  make  it  better  conditioned  while  at  the 
same  time  maintaining  the  singularity  of  those  columns  that  are  used  in  the  calculation  of  the 
determinant.  This  is  accomplished  by  normalizing  all  of  the  columns  of  the  matrix  and 
orthonorm alizing  (using  Gram-Schmidt  orthonormalization)  those  two  that  appear  in  the 
determinant.  If  the  columns  of  the  transfer  matrix  R  are  neither  too  nearly  parallel  nor  too  large 
in  magnitude  it  may  be  possible  to  obtain  a  solution  by  altering  R  in  this  manner.  However,  it 
is  frequently  not  possible  to  accurately  orthonormalize  the  transfer  matrix  for  the  entire  interval; 
in  these  cases  the  interval  of  integration  is  divided  into  two  or  more  smaller  intervals.  The 
solution  pw  at  the  endpoint  of  the  interval  is  then  the  product  of  the  transfer  matrices  for  each  of 
these  intervals  and  the  initial  condition:  for  two  subintervals,  Pw  =  R2  R,  pM.  As  R!  involves 
integration  over  a  shorter  domain,  its  columns  will  not  be  so  large  in  magnitude  nor  so  nearly 
parallel  as  those  of  the  matrix  R;  it  may  thus  be  possible  to  orthonormalize  R„  as  indicated 
above,  to  obtain  a  new  matrix  RL.  This  modified  matrix  is  next  multiplied  by  R2.  This  product 
will  also  have  non-orthonormal  columns,  but  the  resulting  matrix  will  again  be  better  conditioned 
than  the  matrix  R.  It  may  thus  be  possible  to  orthonormalize  the  product  R2  R,  and  use  the 
resulting  matrix  to  evaluate  the  required  determinant.  Clearly  this  procedure  may  be  repeated 
over  as  many  subintervals  as  necessary  to  obtain  accurate  results. 

Note  that  because  c  and  <p  (and  hence  p)  are  complex  valued,  the  Rj  will  likewise  be 
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complex  valued,  and  the  dot  product  used  in  the  orthonormalization  must  therefore  be  applicable 
to  complex  numbers.  Two  possible  choices  are  a*  b  =  a,  b/  (where  the  dagger  indicates 
complex  conjugate,  and  repeated  indices  indicate  summation  over  the  index)  or  a*  b  =  a,  bj.  We 
use  the  second  (complex  valued)  dot  product,  as  it  preserves  the  analyticity  of  the  determinant 
and  allows  the  use  of  a  Muller’s  method  to  locate  the  eigenvalue  (we  use  DZANLY  from  IMSL). 
Some  care  must  be  taken,  however,  as  this  analytic  dot  product  may  be  zero  when  the  vectors 
being  multiplied  are  non-zero,  but  this  never  proved  to  be  a  problem. 

Next  suppose  that  there  are  more  complicated  boundary  conditions  at  the  channel  midline; 
this  occurs,  for  example,  when  we  test  our  code  by  considering  flow  over  a  compliant  plate.  In 
this  case  it  may  not  be  possible  to  select  an  initial  condition  that  corresponds  to  a  linear 
combination  of  basis  vectors,  as  was  the  case  above.  It  will  of  course  still  be  possible  to  write 
down  two  initial  conditions,  but  they  will  in  general  be  linear  combinations  of  all  four  basis 
vectors,  so  that  if  we  continue  as  before  and  orthonormalize  those  columns  of  the  transfer  matrix 
that  contribute  to  the  determinant,  all  four  columns  will  be  orthonormalized.  However,  while  the 
orthonormalization  of  two  columns  of  the  matrix  leaves  the  singularity  of  the  deteiminant 
unchanged,  this  is  not  the  case  for  three  or  more  columns.  To  see  this,  consider  the  determinant 
condition  (5.15),  in  matrix  form 


Det 


Bl'ri  Bi  r4 

B2‘  r2  B2  r4 


=  0. 


(5.16) 


Normalizing  the  column  vectors  ij  will  clearly  not  alter  the  singularity  of  the  determinant. 
Orthogonalization  of  the  second  and  fourth  columns  of  the  transfer  matrix  results  in  the 
replacement  of  r4  with  r4  -  (r2-  rj  r2  /  (r2-  r^;  this  is  equivalent  to  adding  a  multiple  of  the  first 
column  of  the  determinant  matrix  (5.16)  to  the  second,  which  also  leaves  the  singularity  of  the 
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determinant  unchanged.  However,  if  three  columns  of  the  transfer  matrix,  say  r2,  r3  and  r4,  are 
orthogonalized,  then  r3  and  r4  are  transformed  according  to 

r3  -*  r3 A  "  r3  "  Kr3'r2)r2  I  (V  r2 )1  30(3  (5.17) 

r4  -*  r4  -  [(r4-r2)r2  /  (r, •»*,)]  -  [(vrM)rM  /  (rM-rM)] , 

and  orthogonalization  results  both  in  taking  a  linear  combination  of  the  two  columns  of  the 

determinant  matrix  and  in  the  appearance  in  the  determinant  of  factors  of  the  third  column  of  the 

transfer  matrix.  As  the  factors  of  r3  were  not  present  in  the  original  determinant,  we  may  not 

expect  that  the  singularity  of  the  modified  determinant  will  be  that  of  the  original.  This  analysis 

was  simplified  by  the  fact  that  we  left  the  midline  conditions  in  the  form  of  basis  vectors,  but  is 

readily  extended  to  the  more  general  case.  Thus  in  order  to  use  the  orthonormalization  scheme 

developed  above  we  must  be  able  to  begin  with  initial  conditions  in  the  form  of  the  basis  vectors 

ei 

To  accomplish  this,  a  new  variable  q  =  H  p  is  introduced,  where  the  matrix  H  is  chosen 
to  transform  the  boundary  conditions  at  one  end  of  the  interval  to  a  simple  form  such  as 
=  (p'  =  0.  With  the  vector  operators  giving  the  boundary  conditions  at  either  end  of  the  interval 
written  as  Bj  and  Fjt  the  matrix  H  is  thus  defined  so  that  F,*  H'1  =  (1  0  0  0)T  and 
F2-  H1  =  (010  0)T.  This  modified  system  is  then  solved  as  described  above,  bearing  in  mind 
that  the  integration  routine  and  boundary  conditions  Bj  must  be  modified  to  take  into  account  the 
transformation  of  the  variable  p  by  H. 

§  5.4  Comparison  with  known  stability  results 

In  geometries  bounded  by  rigid  plates  we  scale  velocity  on  the  flow  speed  of  the  base 
flow  profile,  and  lengths  on  the  half  channel  width  (for  channels)  or  the  displacement  thickness 
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Comparison  with  previous  stability  results,  TSI 


Case 

Present  work 

Previous  result 

Plane  Poiseuille  Flow 

ReCR  =  5772.2,  k^n  - 
1.021 

ReCR  =  5772.2,  k^R  = 

1.021  1 

Developing  Channel 

Flow,  x=100 

ReCR  =  9383,  k^R  =  1.67 

ReCR  =  9790  2 

x=170 

ReCR  =  8330,  kc„  =  1.34 

ReCR  =  8420  2 

Single  Plate  (Blasius) 

Flow 

ReCR  =  519.7,  kcR  =  0.305 

ReCR  =  520,  kcR  =  0.301  3 

Table  V.i  Comparison  of  critical  Reynolds  and  wave  numbers  for  TSI  in  flows 
with  rigid  boundaries.  Hughes,  in  Drazin  &  Reid  1981  2:  Gupta  &  Garg  1981 
(from  their  figure  1)  3:  Jordinson  1970 

(values  of  for  the  results  of  Gupta  &  Garg  are  difficult  to  obtain  accurately 
from  their  figures  and  so  are  omitted.) 


of  the  boundary  layer  (for  flow  over  a  single  plate).  For  flow  over  a  single  plate  the  channel 
symmetry  or  antisymmetry  conditions  are  replaced  by  the  requirement  that  <p  decay  exponentially 
outside  the  boundary  layer.  The  appropriate  decay  rates  to  require  at  this  point  are  found  from 
consideration  of  the  Orr-Sommerfeld  equation  for  large  z,  in  which  limit  the  equation  has  constant 
coefficients;  solutions  are  thus  decaying  exponentials,  as  found  in  the  exact  solution  for  plug  flow, 
with  decay  rates  =  k  and  m2  =  (k2  +  i  k  (1  -  c)  Re)1/2  (where  rea^m^  >  0  to  avoid  growing 
solutions).  The  boundary  conditions  requiring  that  the  eigenfunction  decay  exponentially  with 


these  decay  rates  are 
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(—  +  k)(J—  -  m£)<p  *  0  and  {—  +  -  *2)<p  =  0  (5.18) 

dz  dz2  dz  dz 2 

(Carpenter  and  Garrad  1985).  As  indicated  above,  these  conditions  require  that  the  problem  be 

reformulated  before  the  orthonormalization  scheme  may  be  applied. 


Ranges  of  unstable  wave  numbers,  unbounded  flow, 
single  compliant  plate 


Reynolds  Number 

present  work,  range  of 
unstable  wavenumbers 

C&G,  range  of  unstable 
wavenumbers 

4000 

0.050  -  0.130 

0.051  -  0.131 

0.142  -  0.344 

0.146  -  0.343 

4500 

0.070  -  0.142 

0.069  -  0.145 

0.158  -  0.380 

0.159  -  0.380 

Table  V.ii  Comparison  of  Blasius  flutter  roots  from  present  and  past  work, 
scaled  on  maximum  flow  speed  and  boundary  layer  displacement  thickness. 
C&G  =  results  of  Carpenter  and  Garrad  (1986,  from  their  figure  13).  E=0.5 
Nmm'2 


In  table  V.I  we  compare  with  previous  work  the  stability  results  we  obtain  for  the  TSI 
for  Poiseuille  and  developing  channel  flow,  and  for  Blasius  flow  over  a  single  plate,  and  in  figure 
5.3  compare  the  neutral  stability  curve  we  obtain  for  Blasius  flow  over  a  single  rigid  plate  with 
that  of  Jordinson  (1970).  The  agreement  is  seen  in  all  cases  to  be  very  good;  the  lack  of  exact 
agreement  with  the  developing  channel  flow  results  of  Gupta  and  Garg  (1981)  is  due  to 
differences  in  the  representation  of  the  developing  flow  profile;  while  we  use  a  perturbation 
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solution  for  the  profile,  they  used  a  finite  difference  solution  of  the  boundary  layer  equations. 

The  numerical  solution  of  the  Orr-Sommerfeld  equation  in  a  compliant  channel  may  be 
checked  by  considering  a  plug  base  flow  and  comparing  with  the  exact  solution  obtained  in 
Chapter  4.  The  exact  and  numerical  solutions  for  this  case  agree  to  the  numerical  precision  used. 

Finally  we  seek  to  duplicate  the  results  of  Carpenter  and  Garrad  (1985, 1986)  for  Blasius 
flow  over  a  compliant  plate.  Their  results  are  also  obtained  using  a  multiple  shooting  algorithm 
with  orthonoimalization,  with  a  polynomial  approximation  being  used  for  the  Blasius  profile.  In 
figure  5.4  points  on  their  neutral  stability  curve  for  the  TSI  (obtained  from  their  figure  11)  are 
shown  with  the  corresponding  curve  from  our  solution;  in  this  figure  we  compare  with  their 
choice  of  wall  elastance  E  =  0.3  Nmm'2  (see  Carpenter  and  Garrad  1985).  We  also  compare  in 
table  V.II  ranges  of  unstable  wave  numbers  for  the  travelling  wave  (Type  B)  flutter  instability 
(Carpenter  and  Garrad  1986).  For  either  instability  the  agreement  between  the  present  work  and 
previous  results  is  very  good. 

§5.5  Results 

The  flutter  dispersion  relation  found  numerically  for  developing  flow  is  shown  as  a 
function  of  the  flow  speed  S  in  figure  5.5.  As  with  the  plug  flow  base  state  considered  in 
Chapter  4,  there  are  roots  travelling  downstream  and  upstream  for  small  S,  which  we  call  roots 
1  and  2  respectively;  root  1  is  shown  in  figure  5.5.  Root  2,  as  with  the  plug  flow,  slows  and 
reverses  direction  as  S  is  increased,  becoming  unstable  only  after  the  appearance  of  instability  for 
root  1.  The  dispersion  relation  in  5.5  is  for  symmetric  disturbances;  that  for  antisymmetric 
disturbances  is  qualitatively  similar.  Neutral  stability  curves  for  both  symmetric  and 
antisymmetric  disturbances  are  shown  in  figure  5.6,  showing  the  symmetric  mode  to  be  the  least 
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Figure  5.3:  Neutral  stability  curve  for  Blasius  flow  over  a  single  rigid  plate,  TSI. 
o  present  work;  •  Jordinson  (1970) 
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Figure  5.4:  Neutral  stability  curve  for  Blasius  flow  over  a  single  compliant  plate, 
TSI.  Dashed  curve,  present  work;  •  Carpenter  and  Garrad  (1985). 
(Elastance  =  0.3Nmm'2) 
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Figure  5.6:  Neutral  stability  curves,  flutter  instability.  Solid  curve,  symmetric 
mode;  Dashed,  antisymmetric.  Solid  dot  gives  long  wave  instability 
(Scrlw  );  open  dot,  flutter  (SCRF ).  M  =  3190,  B  =  542,  T  =  0,  G  =  0.054, 
d  =  882,  =  2230,  Xo  =  2.0. 
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stable.  The  long  wave  limit  is  shown  on  the  symmetric  curve  as  a  solid  dot;  as  shown  in  §5.2, 
there  is  no  long  wave  instability  for  antisymmetric  disturbances.  We  denote  the  absolute 
minimum  of  the  (symmetric)  neutral  stability  curve  as  SCR  and  the  corresponding  value  of  the 
wave  number  as  kcR.  The  critical  flow  speed  for  the  long  wave  instability  (solid  dot)  we  identify 
as  Scrlw,  and  the  local  minimum  for  k*0  (shown  in  5.6  with  an  open  dot)  as  SCRF,  the  critical 
flow  speed  for  flutter.  In  figure  5.6  SCR  =  SCRF.  Comparison  of  the  developing  flow  neutral 
stability  curves  in  figure  5.6  with  those  for  plug  flow  in  figure  4.6  of  Chapter  4  shows  the 
qualitative  difference  that  the  long  wave  instability  present  in  the  developing  flow  is  absent  in  the 
plug  flow.  The  developing  flow  is  also  significantly  less  stable  than  the  plug  flow. 

As  was  noted  in  §5.2,  because  the  critical  flow  speed  for  the  long  wave  instability 
depends  on  wall  characteristics  such  as  the  mass  ratio  M,  etc.,  only  at  higher  orders  in  k,  it  is 
possible  by  varying  system  parameters  to  stabilize  the  finite  wavelength  instability  and  make  the 
long  wave  instability  critical.  This  is  shown  in  figure  5.7,  in  which  neutral  stability  curves  are 
shown  for  a  number  of  different  mass  ratios.  Reduction  of  the  mass  ratio  M  (e.g.  by  reducing 
the  wall  density)  stabilizes  the  flutter  instability  (increases  SCRF ),  so  that  the  instability  moves 
from  non-zero  to  zero  wave  number  as  M  is  decreased.  Once  the  long  wave  instability  has 
appeared,  further  variation  in  M  does  not  alter  the  onset  flow  speed  of  the  instability,  as  discussed 
in  §5.2  and  below.  This  transition  between  the  finite  and  long  wavelength  instabilities  admits  the 
appearance  of  a  co-dimension  two  bifurcation  point,  occurring  when  the  long  wave  and  flutter 
instabilities  have  the  same  critical  flow  speed  (SCR  =  SCRLW  =  SCRF). 

In  figure  5.8A  the  effect  of  varying  the  wall  damping  G  and  mass  ratio  M  on  the  system 
is  shown  by  plotting  SCRLW  and  SCRF  as  functions  of  G  for  two  values  of  M.  As  G  is  increased 
or  M  decreased,  the  system  is  stabilized  with  respect  to  flutter  (SCRF  increases),  but  the  long  wave 
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Figure  5.7:  Neutral  stability  curves  for  different  mass  ratios,  showing  long  wave 
and  flutter  instabilities.  Solid  dot  gives  long  wave  limit,  o  M  =  500; 
v  M  =  100;  *  M  =  10;  □  M  =  1.  B  =  542,  T  =  0,  Gx5000  =  0.0845, 
d  =  882, =  2230,  Xo  =  1.0.  (6  ■  b*  G*  /  fl) 
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Figure  5.8:  Critical  flow  speeds  SCRLW  and  SCRF  (5.8A)  and  frequency  at  SCR 
(5.8B)  as  functions  of  wall  damping  and  mass  ratio.  Dotted  line  gives 
Scrlw;  Symbols  give  SCRF  in  5. 8 A  and  frequency  in  5.8B:  o  M  =  5000; 
•  M  =  500.  B  =  542,  T  =  0,  d  =  882, =  2230,  Xq  =  1.0. 
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limit  (Scrlw)  is  unchanged,  so  that  for  large  wall  damping  stability  is  lost  to  the  long  wave 
instability.  The  stabilization  of  the  flutter  instability  is  consistant  with  its  being  Type  B,  and 
being  stabilized  by  reduction  of  wall  inertia.  Note  that  5.8A  is  qualitatively  different  than  the 
behavior  of  the  plug  base  flow  shown  in  figure  4.7A  and  discussed  in  §4.4,  for  which  the  effect 
of  increasing  M  could  be  stabilizing  for  large  G  because  of  the  plug  flow  long  wave  stability. 
Figure  5.8B  shows  the  frequency  of  the  oscillations  appearing  for  developing  flow  at  SCR.  The 
behavior  of  the  frequency  is  similar  to  that  seen  for  plug  flow  (figure  4.7B),  decreasing  with 
increasing  M  or  G,  until  the  appearance  of  the  long  wave  instability  results  in  a  discontinuous 
drop  to  zero  frequency. 

Figure  5.9A  shows  SCRF  versus  wall  elastance  E  (nondimensionalized  on  (p*v*2/ 
b‘3)xia*)  for  two  values  of  the  wall  flexural  rigidity  B.  For  the  parameter  values  considered, 
SCR=  Scrf.  The  frequency  of  the  oscillations  appearing  at  SCR  is  shown  in  figure  5.9B. 
Stiffening  the  wall  by  increasing  either  E  or  B  increases  both  SCR  and  the  corresponding  frequency 
of  oscillations;  this  is  similar  to  the  results  for  plug  flow  (figure  4.8)  and  is  expected  physically. 
Note  that  the  frequency  in  5.9B,  as  that  in  4.8B,  increases  approximately  as  E1/2,  consistent  with 
the  results  of  Grotberg  and  Reiss  (1984).  The  inset  plots  in  5.9  show  the  difference  between  the 
values  of  SCRF  for  the  two  values  of  B  (in  5.9A)  and  the  corresponding  difference  in  frequency 
(in  5.9B),  showing  that  the  effect  of  B  is  larger  for  larger  wall  elastance,  for  which  the  instability 
has  shorter  wavelength.  This  is  as  expected,  and  a  similar  effect  may  be  obtained  through 
variation  of  the  wall  tension  T.  Note  that  the  second  wall  stiffness  parameter,  d,  is  not 
independent  of  B,  varying  only  by  a  factor  of  12  (b*  /  h*)2;  we  have  retained  it  as  a  second 
parameter  only  to  facilitate  comparison  with  the  results  obtained  for  the  limit  of  no  horizontal 
wall  motion  (which  corresponds  to  d'1  —>  0).  We  find  little  difference  between  the  cases  d'1  =  0 
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Figure  5.9:  Critical  flow  speed  SCR  (=  SCRF)  (5 .9 A)  and  flutter  frequency  (5.9B) 
as  functions  of  wall  elastance  E.  Elastance  scaled  on  (p*  v*2  /  b‘3 )  x  104, 
velocity  on  (v*  /  b*  )xl02,  frequency  on  (v*  /  b* 2 ).  o  B  =  2.688  x  104; 
•  B  =  2.688  x  10s.  Inset  plots  give  the  difference  between  SCR  and  fcR 
for  the  two  values  of  B,  as  a  function  of  E.  Other  parameters  as  figure 
5.7  (M  =  3190),  with  modified  scaling. 
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and  d'1  *  0. 

Because  the  half-channel  width  b*  also  appears  in  the  scaling  of  the  problem,  to  examine 
the  effect  of  changing  b*  we  scale  the  system  as  noted  in  §5.2,  scaling  lengths  on  and 
velocities  on  (v*  /  Xq*  )xl02.  Figure  5.10  shows  the  neutral  stability  curves  for  a  number  of 
different  half  channel  widths  b;  as  b  is  increased,  the  effect  of  the  second  plate  and  hence  of 
viscosity  in  the  channel  is  decreased,  which  stabilizes  the  Type  B  flutter.  However,  as  indicated 
in  §5.2,  the  long  wave  instability  continues  to  appear  at  lower  flow  speeds,  so  that  for  sufficiently 
wide  channels  instability  will  be  lost  to  the  long  wave  instability.  This  is  seen  in  figure  5.1, 
which  shows  SCRLW  and  SCRF  as  functions  of  b;  we  see  that  as  b  is  increased  past  =9  the  instability 
goes  from  the  flutter  to  the  long  wavelength  instability.  This  is  particularly  interesting  in  light 
of  the  experimental  observations  (Gavriely  et.al.  1989)  that  flutter  only  occurs  following  tube 
collapse;  we  hypothesize  that  this  tube  collapse  is  the  physical  manifestation  of  the  long  wave 
instability  and  that  flutter  is  then  seen  only  when  the  channel  width  is  reduced  enough  that  it 
may  become  critical.  As  the  wall  elastance  (Rw)  is  increased,  the  codimension  two  point  moves 
to  larger  channel  widths.  This  is  consistent  with  physical  observations  of  flutter  in  stiffer  tubes, 
which  do  not  collapse;  the  shift  in  the  location  of  the  codimension  two  point  in  this  case  is 
evidently  sufficient  to  ensure  the  appearance  of  the  flutter  instability.  With  increasing  half¬ 
channel  width,  the  perturbation  solution  for  the  developing  flow  that  we  are  using  converges  to 
Blasius  flow  over  a  single  compliant  plate,  and  the  stability  results  for  the  single  plate  system 
(e.g.  Carpenter  and  Garrad  1986)  are  recovered  for  the  finite  wave  number  instability.  However, 
for  long  waves  the  effect  of  the  opposing  channel  wall  is  always  significant,  so  that  the  single 
plate  limit  is  not  reached;  this  is  shown  in  figure  5.10  through  the  inclusion  of  a  stability  curve 
for  Blasius  flow  over  a  single  compliant  plate.  We  discuss  the  effect  of  increasing  channel  width 


80 


3 

Flow 

Speed 

S 


2 


1 

0.00  0.01  0.02  0.03  0.04 

Wave  Number  k 


Figure  5.10:  Neutral  stability  curves  for  different  channel  widths,  showing  long 
wave  instability  and  comparison  with  Blasius  flow  over  a  single 
compliant  plate.  Lengths  scaled  on  Xo‘,  velocities  on  (v*  /  x0*)xl02. 
o  b  =  4.82;  •  b  =  14.5;  v  b  =  24.1;  ▼  b  =  72.4;  dotted  curve  is  for 
Blasius  flow  over  a  single  compliant  plate.  Other  parameters  as  figure 
5.7  (M  =  3190),  with  modified  scaling. 
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in  slightly  greater  detail  in  §5.6.  Figure  5.1  A  may  also  be  compared  with  figure  4.9A  of  Chapter 
4,  which  shows  SCRF  for  plug  flow  as  a  function  of  half  channel  width.  The  behavior  of  the 
critical  flow  speed  is  similar  for  either  base  flow,  except  for  the  long  wave  instability  which  is 
absent  for  plug  flow.  Comparing  5.1B  and  4.9B,  the  change  in  flutter  frequency  is  also  similar 
for  narrower  channels;  however,  for  wider  channels,  the  effect  of  the  long  wave  instability  is  seen 
for  the  developing  flow  while  the  plug  flow  long  wave  stability  increases  the  plug  flow  frequency 
as  discussed  in  §4.4. 

Because  increasing  the  Reynolds  number  of  the  flow  is  equivalent  to  decreasing  the  fluid 
viscosity,  such  an  increase  stabilizes  the  flutter  instability.  However,  in  the  case  of  the  developing 
profile,  increasing  the  Reynolds  number  also  changes  the  base  flow  profile,  causing  (for  constant 
flow  speed  S  and  axial  position  Xq  )  a  narrowing  of  the  boundary  layer  at  the  wall.  A  similar 
effect  may  be  obtained  by  decreasing  Xq  ;  this  destabilizes  the  flutter  instability.  It  is  difficult  to 
draw  conclusions  for  wide  ranges  of  axial  positions,  however,  due  to  limitations  imposed  by  the 
downstream  limit  of  validity  for  our  developing  profile.  Further,  as  points  farther  upstream  are 
considered,  non-parallelism  in  the  flow  may  become  significant;  we  return  to  this  issue  in  the 
discussion  in  Chapter  7. 

As  noted  in  Chapter  4,  the  effects  that  we  find  destabilize  the  flutter  instability  correspond 
with  the  distinguishing  characteristics  of  the  lung  airways  of  wheezing  patients.  For  developing 
flow,  however,  there  is  also  the  long  wave  instability,  which  we  expect  is  related  to  tube  collapse, 
and  hence  flow  limitation.  As  wide  channels  are  more  susceptible  to  the  long  wave  instability 
than  to  flutter,  this  agrees  with  clinical  observations  that  forced  expiratory  wheezes  are  only  seen 
following  flow  limitation  (Gavriely  et.al.  1987).  In  Chapter  6  we  compare  the  theoretical  results 
from  the  developing  flow  model  with  experimental  results. 
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§5.6  Discussion  of  the  limit  of  infinite  channel  width 

To  consider  results  for  very  wide  channels,  it  is  convenient  to  investigate  the  effect  of 
increasing  channel  width  on  the  flutter  dispersion  relation  shown  in  figure  5.5.  In  this  figure  the 
eigenvalue  c  is  seen  to  change  smoothly  with  increasing  flow  speed;  however,  as  the  half  channel 
width  b  is  increased,  there  appears  near  the  point  at  which  the  base  flow  speed  becomes  equal  to 
the  phase  speed  of  the  disturbance  an  abrupt  drop  in  the  phase  speed  and  upward  bend  in  the 
growth  rate  of  the  disturbance.  As  b  is  increased  still  further,  this  becomes  an  actual 
discontinuity;  this  is  shown  in  figure  5.1 1A.  To  illustrate  this  more  clearly,  figure  5.1  IB  shows 
the  eigenvalue  in  the  complex  c  plane;  in  this  form,  imag(c)  is  given  as  a  function  of  real(c),  with 
flow  speed  S  parameterizing  the  curves  (only  the  region  of  parameter  space  about  S  =  real(c)  is 
shown).  The  dependence  of  this  behavior  on  the  condition  S  =  real(c)  is  emphasized  by  scaling 
c  on  the  maximum  flow  speed  of  the  base  flow,  so  that  the  line  S  =  real(c)  lies  along  real(c)  = 
1.  Similar  behavior  may  be  obtained  by  varying  other  system  parameters;  this  occurs  when  the 
effective  channel  width  is  very  large  {e.g.  the  mass  ratio  M  goes  as  (b  )  ',  so  that  decreasing  M 
may  be  interpreted  as  increasing  the  channel  width,  subject  to  variation  of  other  parameter 
values). 

To  understand  this  behavior,  it  useful  to  consider  more  specifically  the  characteristics  of 
the  Orr-Sommerfeld  system  in  bounded  and  unbounded  domains.  In  an  unbounded  domain, 
unlike  the  finite  case,  there  are  only  a  finite  number  of  eigenvalues  in  the  discrete  spectrum 
(MiklavCiC  1983),  and  these  are  supplemented  by  a  continuous  spectrum  consisting  of  ’improper’ 
eigenvalues  having  phase  speeds  equal  to  the  maximum  flow  speed  of  the  base  profile  and 
negative  growth  rates  (Grosch  and  Salwan  1978;  Craik  1991).  As  the  Reynolds  number  of  the 
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flow  is  increased,  eigenvalues  are  added  to  the  discrete  spectrum,  appearing  sequentially  at  the 
line  in  the  complex  c  plane  giving  the  continuous  spectrum  (Mack  1976).  For  Blasius  flow  over 
a  rigid  plate,  the  only  eigenvalues  in  the  vicinity  of  the  continuous  spectrum  are  significantly 
damped,  but  this  is  not  the  case  when  a  flexible  wall  is  considered.  In  this  case  the  Type  B 
flutter  mode  (which  is  our  root  1)  may  be  critical,  and  may  have  a  phase  speed  that  approaches 
the  speed  of  the  base  flow.  The  presence  of  the  continuous  spectrum  then  complicates  matters, 
as  it  is  difficult  to  accurately  locate  eigenvalues  near  the  continuous  spectrum  (Carpenter  and 
Garrad  1986).  In  a  finite  domain  there  is  no  continuous  spectrum  (Lin  1961),  and  this  difficulty 
does  not  arise,  but  as  the  finite  domain  is  modified  to  look  more  like  the  infinite  case  (when  the 
channel  width  is  increased,  or  other  parameters  varied  to  obtain  a  similar  result)  it  may  be 
expected  that  the  numerical  solution  to  the  problem  might  encounter  difficulties  as  the  phase 
speed  of  the  disturbance  becomes  similar  to  the  flow  speed  of  the  base  flow,  as  seen  in  figure 
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CHAPTER  6:  COMPARISONS  WITH  EXPERIMENTS 


§6.1  Overview 

As  the  motivation  for  this  study  is  to  model  the  oscillatory  instabilities  in  the  lung  airways 
that  lead  to  wheezing  lung  sounds,  we  would  like  to  compare  our  theoretical  predictions  with 
actual  experimental  observations  of  such  wheezing.  However,  the  characteristics  of  the  lung 
airways  that  we  need  to  obtain  the  parameters  in  our  model  (M,  B,  etc.)  are  not  known  with  any 
accuracy,  making  such  comparison  difficult.  We  therefore  first  consider  the  theoretical  modelling 
of  experiments  with  flexible  tubes,  in  §6.2  -  §6.4,  and  return  to  the  issue  of  the  lung  in  §6.5.  In 
that  Gavriely  et.al.  (1987,  1989)  have  demonstrated  that  the  characteristics  of  both  the  collapse 
phenomenon  and  the  sounds  produced  in  such  tubes  are  similar  to  those  in  the  lung,  this 
comparison  may  itself  indicate  that  our  model  is  applicable  to  the  lung  and  shed  insight  on  the 
lung  experiments. 

§6.2  Comparison  with  flutter  in  tube  experiments 

Gavriely  et.al.  (1989)  experimentally  investigated  air-conducting  flexible  tubes,  and  found 
an  oscillatory  instability  that  followed  tube  collapse.  Based  on  the  characteristics  they  give  for 
their  tubes  and  flow  (h‘=0.19cm,  b*=0.07cm  (after  collapse,  at  the  onset  of  flutter), 
pw*=0.941g/cm3,  p‘=8xl04g/cm\  v‘=0.225cm2/s,  T*=0,  G’=  1.72s1,  E'=5. 86x10s  dyne/cm3,  and 
D*=Y*h‘3/9  (where  Y*=the  Young’s  modulus)=7621dyne  cm),  we  find  the  following  non- 
dimensional  parameter  values  for  our  model:  M  =  3190,  B  =  542,  T  =  0,  G  =  0.0542,  d  =  882, 
and  =  2230.  For  the  experimental  case  for  which  these  parameter  values  were  derived,  the 
critical  flow  speed  at  which  the  oscillatory  instability  was  first  observed  was  2670cm/s;  this  value 
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was  obtained  by  averaging  the  flow  rate  over  the  cross-sectional  area  of  the  tube.  The  symmetric 
instability  appearing  at  this  flow  speed  had  a  frequency  of  300Hz. 

To  obtain  results  using  our  developing  flow  model,  it  is  necessary  not  only  to  use  the 
parameter  values  indicated  above,  but  also  to  select  an  axial  position  at  which  to  evaluate  the 
developing  flow  profile.  In  the  experiments  a  long  (60cm)  flexible  tube  was  fastened  at  either 
end  to  a  rigid  pipe,  and  air  flow  through  the  tube  was  induced  by  a  suction  pump  at  the 
downstream  end.  As  the  driving  pressure  gradient  was  increased,  collapse  and  flutter  occurred 
at  the  downstream  end  of  the  tube.  The  development  length  for  a  tube  is 
xent  =  0.25  r* 2  uavg*  /  v*  (where  r*  is  the  tube  radius  and  uAVG*  the  average  fluid  velocity) 
(Schlichting  1955);  for  the  experiment  the  uncollapsed  tube  radius  r*  =  0.325cm,  v*  =  0.225cm2/s, 
and  uavg*  =  1600cm/s,  so  that  x^  =  188cm.  Thus  it  is  clearly  appropriate  to  consider  the  flow 
in  the  collapsed  tube  section  to  be  developing.  However,  the  precise  appearance  of  the  profile 
is  unclear,  as  when  the  flow  enters  the  collapsed  section  it  regains  a  more  pronounced  boundary 
layer  character  due  to  the  constriction.  We  therefore  calculate  the  stability  of  the  system  at 
several  axial  locations,  and  also  compare  with  the  results  for  fully  developed  flow;  an  additional 
comparison  is  then  made  with  the  predictions  of  the  plug  flow  solution. 

The  predicted  flutter  frequencies  for  a  developing  flow  at  axial  positions  x  =  1,  2  and  3 
half  channel  widths  (recall  that  it  is  necessary  that  we  consider  Xo  small  for  the  for  the  developing 
profile  to  be  valid;  this  requirement  is  slightly  less  stringent  than  x^  <  3)  are  312Hz,  31 1Hz,  and 
310Hz;  for  Poiseuille  flow,  308Hz;  and  for  plug  flow,  333Hz.  The  numerical  predictions  for  the 
viscous  (developing  or  Poiseuille)  profiles  are  thus  all  quite  close  to  the  experimental  value,  while 
the  plug  flow  result  is  higher.  That  is  should  be  larger  than  the  numerical  values  is,  however, 
reasonable  from  the  trend,  seen  in  the  frequencies  for  the  developing  flow,  of  increasing 
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frequency  with  decreasing  axial  position. 

The  critical  flow  speeds  SCRF  corresponding  to  the  frequencies  above  are  for  the 
developing  flow  344cm/s,  417cm/s  and  463cm/s  (at  x<,  =  1,  2,  and  3);  for  Poiseuille  flow, 
540cm/s;  and  for  plug  flow,  l,488cm/s.  It  is  not  appropriate  to  directly  compare  the  theoretical 
and  experimental  flow  speeds,  however,  as  the  value  measured  experimentally  (2670cm/s)  is  an 
average  over  a  collapsed  tube  cross-section.  As  this  cross-section  is  a  dumbbell-like  shape,  the 
velocity  in  the  central  (channel)  section  of  the  collapsed  tube  (where  the  oscillations  appear)  will 
be  lower  than  that  in  the  outer  lobes  because  of  viscous  resistance. 

To  estimate  the  degree  to  which  the  channel  flow  speed  might  differ  from  the 
experimentally  measured  cross-sectional  average,  we  use  the  finite  elements  package  FIDAP  to 
calculate  the  velocity  profile  in  a  collapsed  tube  cross-section.  This  requires,  however,  knowledge 
of  the  dimensions  of  the  collapsed  tube  section,  which  are  not  reported  by  Gavriely  et.al  (1989); 
dimensions  that  are  given  are  the  half-channel  width  (0.07cm)  and  cross-sectional  area  (0.2cm2 ). 
Flaherty  et.al.  (1972)  numerically  obtained  shapes  for  collapsed  tubes,  to  which  the  experimental 
dimensions  might  be  matched,  but  their  model  is  for  a  tube  with  walls  of  negligible  width  while 
in  the  experiments  the  width  of  the  tube  wall  (0.19cm)  was  a  substantial  percentage  of  the 
uncollapsed  radius  (0.325cm).  We  therefore  consider  two  approximate  collapsed  tube  shapes 
based  on  the  dimensions  that  are  given  and  the  results  of  Flaherty  et.al.',  the  two  differ  in  that  the 
second  is  flatter,  with  smaller  side  lobes  and  a  longer  channel  section.  By  considering  the  two, 
we  obtain  a  qualitative  estimate  of  how  the  tube  shape  may  influence  the  flow  speeds  in  the  cross- 
section.  As  the  imprecision  in  estimating  the  tube  shape  and  dimensions  itself  prevents  the 
velocity  calculation  from  being  quantitative,  only  fully  developed  flow  is  considered,  as  this 
results  in  an  order  of  magnitude  decrease  in  the  required  computational  resources  from  those 
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necessary  for  the  calculation  of  developing  flow.  Details  of  the  FIDAP  calculation  appear  in 
Appendix  B.  For  the  flow  rate  of  the  experiment  (533cm3/s),  the  average  velocity  in  the  channel 
section  of  the  two  shapes  is  found  to  be  1440cm/s  and  lOlOcm/s,  respectively.  These  calculations 
thus  indicate  that  the  velocity  at  the  center  of  the  collapsed  tube  section  may  be  significantly 
lower  than  the  cross-sectionally  averaged  value  determined  in  the  experiments,  and  thus  that  the 
actual  channel  flow  speed  in  the  experiments  was  likely  to  have  been  closer  to  our  theoretical 
predictions. 

An  additional  complication  in  this  comparison  arises  from  the  geometry  of  the  collapse 
itself,  which  results  in  both  a  constriction  (in  the  transverse  dimension  of  the  channel)  and 
expansion  (out  of  the  channel,  into  the  side  lobes  of  the  dumbbell)  of  the  tube.  The  TSI  is 
stabilized  for  flow  into  a  constriction  (White  1974),  but  the  effect  of  a  lateral  expansion  (and  the 
resulting  three  dimensionality  of  the  flow)  has  not  to  our  knowledge  been  investigated  even  for 
the  TSI  Get  alone  flutter).  Neither  of  these  effects  are  included  in  our  model,  and  they  may, 
along  with  the  variation  in  flow  speeds  investigated  above,  reconcile  the  difference  between  the 
theoretical  and  experimental  critical  flow  speeds. 

§6.3  Comparison  with  collapse  in  tube  experiments 

Next,  we  seek  to  compare  the  experimentally  observed  flow  speed  at  which  tube  collapse 
began  with  the  long  wave  onset  flow  speed  we  find  theoretically.  It  is  difficult  to  determine  the 
experimental  value  from  Gavriely  et.al.  (1989),  as  they  did  not  report  the  critical  flow  speeds  for 
tube  collapse.  However,  we  estimate  the  beginning  of  collapse  to  be  when  the  pressure  -  flow- 
rate  relationship  they  show  first  deviates  from  the  slow  linear  increase  seen  before  collapse  and 
flow  limitation.  From  their  figure  3.B,  this  gives  a  critical  flow  rate  of  approximately  27  1/min, 
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or  an  average  flow  speed  SCRLW"P  (recalling  that  the  uncollapsed  tube  radius  was  0.325cm)  of 
1340cm/s.  Reading  from  our  figure  5.1,  we  find  for  this  channel  width  that  the  long  wave 
instability  is  critical,  with  SCRLW  of  about  680cm/s,  or  roughly  half  of  the  estimate  from  the 
experiments.  Thus  the  developing  flow  model  accurately  predicts  the  existence  of  the  long  wave 
(collapse)  instability.  In  that  the  comparison  between  the  theoretical  and  experimental  critical 
flow  rates  is  at  best  qualitative,  due  to  the  nature  of  our  experimental  estimate  and  the  fact  that 
before  collapse  it  may  be  less  appropriate  to  model  the  tube  as  a  channel,  we  conclude  that  there 
is  evidence  to  support  our  hypothesis  that  tube  collapse  is  related  to  the  long  wave  instability. 

§6.4  The  TSI  in  the  tube  experiments 

Two  observations  verify  that  the  flutter  instability  investigated  above  is  that  which  is 
relevant  to  the  tube  experiments.  First,  the  experimentally  observed  instability  is  symmetric, 
while  the  TSI  is  antisymmetric,  and  second,  we  also  numerically  calculate  the  critical  flow  speed 
for  the  TSI  to  demonstrate  that  it  occurs  for  much  higher  flow  speeds.  This  critical  Reynolds 
number  (scaled  on  maximum  flow  velocity  and  half  channel  width)  is  ReCR  =  5761,  which  is 
equivalent  to  an  average  flow  speed  of  approximately  18,500cm/s,  which  is  far  in  excess  of  both 
the  experimentally  measured  and  theoretically  predicted  critical  flow  speeds  for  flutter.  The 
frequency  of  the  oscillation  for  the  instability  is  11,360Hz  at  criticality,  which  is  again  several 
orders  of  magnitude  larger  than  the  experimental  or  theoretical  values.  These  calculations  are  for 
Poiseuille  flow;  the  critical  Reynolds  numbers  for  the  TSI  in  developing  flows  are  much  higher 
(Chen  and  Sparrow  1967). 
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§6.5  Comparison  with  wheezing  in  the  lung 

We  now  turn  to  the  issue  of  modelling  wheezing  in  the  lung.  For  this  comparison  we 
consider  experiments  with  forced  expiratory  wheezes  (Gavriely  et.al.  1987)  for  which  wheeze 
frequencies  and  flow  rates  were  reported.  The  subjects  in  this  study  executed  forced  expiratory 
maneuvers  to  produce  wheezing  sounds;  a  controlled  suction  pump  was  used  to  facilitate  the 
production  of  the  wheezes,  which  had  frequencies  of  580-2730Hz.  The  volume  flow  rate  of  the 
expirations  was  on  the  order  of  1 1/s.  Wheezes  were  measured  with  a  tracheal  microphone,  but 
the  actual  location  of  the  wheezes  in  the  lung  is  not  indicated.  Kramen  (1983)  also  investigated 
forced  expiratory  wheezes,  though  without  mechanical  forcing,  and  found  them  to  occur  between 
the  trachea  and  lobar  bronchi.  The  frequencies  of  the  wheezes  reported  by  Kramen  are  similar 
to  those  of  Gavriely  et.al.,  as  are  those  in  the  limited  forced  expiratory  experiment  of  Hinchey 
and  Snellen  (1987). 

To  apply  our  theoretical  models  to  these  experiments,  we  need  estimates  for  the  physical 
characteristics  of  the  airways  in  question.  Dragon  and  Grotberg  (1991)  suggest  for  the  trachea 
r*=0.9cm  (uncollapsed  radius),  h‘=0.05cm,  pw*=0.997g/cm3,  G*=5s1,  and  E*  up  to  2.8xl07 
dyne/cm2,  with  v*=0.15cm2/s  and  p*=0.001 14g/cm3.  For  the  bronchi,  we  expect  the  elastance  E* 
to  be  lower,  and  the  radii  r*  to  be  smaller,  Olson  et.al.  (1970)  give  r*=0.65cm,  0.47cm  and  0.36cm 
for  the  first  three  generations  of  bronchi.  To  estimate  a  Youngs  modulus  Y*  for  the  trachea,  and 
collapsed  ’channel-width’  for  the  airways,  we  assume  that  the  ratios  Y*/E*  (=200)  and  b7r* 
(=0.215)  for  tube  experiments  (Gavriely  et.al.  1989)  may  be  applied  to  the  lung.  We  find  E*  and 
Y*  for  the  bronchi  by  decreasing  the  corresponding  values  for  the  trachea  by  a  factor  of  5  for  each 
airway  generation.  Values  of  h*  for  the  bronchi  are  found  by  letting  the  ratio  of  h*  to  r*  be  the 
same  as  for  the  trachea.  Finally,  we  take  for  the  trachea  E*=lxl06  dyne/cm2.  Our  estimates  for 
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parameter 

Trachea 

1°  Bronchi 

2°  Bronchi 

3°  Bronchi 

* 

r 

0.90cm 

0.65cm 

0.47cm 

0.36cm 

b* 

0.19cm 

0.14cm 

0.10cm 

0.077cm 

h* 

0.05cm 

0.04cm 

0.03cm 

0.02cm 

(1)  Pw* 

0.997g/cm3 

0.9g/cm3 

0.9g/cm3 

0.9g/cm3 

(2)  G* 

5s'1 

5s’1 

5s-1 

5s-1 

E*  (dyne/cm2) 

lxlO6 

2x10s 

4X104 

8xl03 

Y*  (dyne/cm2) 

2x10* 

4xl07 

8xl06 

1.6X106 

a 

29,620cm/s 

13,250cm/s 

5,923cm/s 

2,649cm/s 

M 

230 

226 

237 

205 

B 

2.13 

3.70 

6.00 

5.06 

(3)  T 

0 

0 

0 

0 

G 

7.38x1  a3 

0.0119 

0.0200 

0.0298 

d 

369 

544 

800 

899 

Rw 

37,515 

12,362 

3,950 

1,360 

Table  VI.i:  Dimensional  and  non-dimensional  parameters  for  lung  airways.  (1): 
density  of  the  bronchi  is  assumed  to  be  somewhat  less  than  that  of  the  trachea, 
but  of  similar  magnitude.  (2):  wall  damping  is  arbitrarily  taken  to  be  constant 
for  all  airway  generations.  (3):  we  initially  assume  zero  longitudinal  wall  tension 
for  all  airways. 


the  airway  characteristics  and  the  corresponding  dimensional  and  non-dimensional  parameters 
used  in  our  models  are  summarized  in  table  VI.i. 

In  table  Vl.ii  we  present  results  from  the  developing  and  plug  flow  models  of  Chapters 
4  and  5  for  the  parameter  values  indicated  in  table  VI.i.  Because  of  the  crude  nature  of  some  of 
the  estimates  leading  to  those  parameter  values,  we  also  show  results  for  some  other  values.  The 
values  for  the  flutter  frequencies  in  the  table  may  be  compared  with  the  experimental  values  (580- 
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parameter  set: 

Trachea 

1°  Bronchi 

2°  Bronchi 

3°  Bronchi 

f  * 

rCR 

Devel. 

Flow  Scrf  * 

0*0=2) 

o  LW* 
°CR 

1646Hz 

1011Hz 

619Hz 

386Hz 

2308cm/s 

1183cm/s 

593cm/s 

304cm/s 

2395cm/s 

1634cm/s 

878cm/s 

530cm/s 

Plug  fCR* 

Flow 

s  F* 

&CR  I 

1808Hz 

1117Hz 

687Hz 

434Hz 

9544cm/s 

4350cm/s 

1934cm/s 

845cm/s 

parameter  set: 

2°  Bronchi1 

2°  Bronchi* 

2°  Bronchi*1 

2°  Bronchi** 

f  * 

*CR 

Devel. 

Flow  Scrf  * 

0*o=2) 

o  LW* 
°CR 

672Hz 

896Hz 

508Hz 

621Hz 

399cm/s 

568cm/s 

547cm/s 

420cm/s 

878cm/s 

878cm/s 

609cm/s 

609cm/s 

Plug  fCR* 

Flow 

o  F* 
^>CR 

765Hz 

1116Hz 

607Hz 

763Hz 

1236cm/s 

2005 cm/s 

2204cm/s 

1517cm/s 

Table  VI. ii:  Theoretically  predicted  flutter  frequencies  and  critical  flow  rates  for 
the  parameters  indicated  in  table  VLi.  t.$:  parameter  values  as  2°  Bronchi,  table 
VI. i,  with  t:  G  =ls  ‘;  $:  G*=ls1,  T=10;  tt:  E*=2.041xl04,  T=10; 

E*=2.041xl04,  T=10. 


2730Hz),  and  are  seen  to  agree  quite  well.  As  noted  in  §6.2,  it  is  difficult  to  compare  our 
theoretically  predicted  critical  flow  speeds  with  experimental  values,  and  this  difficulty  is  even 
greater  for  the  lung  experiments,  for  which  even  the  collapsed  tube  cross-sectional  area  is  not 
known.  It  is  worth  noting,  however,  that  flutter  (rather  than  the  long  wave  collapse)  is  predicted 
for  the  developing  flow,  consistent  with  our  modelling  already  collapsed  airways.  If  we  assume 
that  the  ratio  of  collapsed  to  uncollapsed  tube  cross-sectional  area  (AcOL  /  7tr* 2 )  is  similar  for  the 
lung  airways  and  tube  experiments  (of  Gavriely  et.al.  1989),  we  have  Aco^l^cm2  for  the 
trachea,  and  =0.80cm2, 0.42cm2  and  0.41cm2  for  the  first  three  generations  of  bronchi.  For  all/s 
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expiratory  flow  rate,  then,  the  average  flow  speed  through  a  collapsed  trachea  would  be  653cm/s; 
through  the  bronchi  (noting  that  there  are  2, 4,  and  7  bronchi  in  the  first  3  generations  of  bronchi 
(Olson  et.al.  1970)),  625cm/s,  595cm/s  and  348cm/s,  respectively.  Based  on  the  comparison  with 
flutter  in  tube  experiments  (§6.2)  and  the  imprecision  in  these  estimates,  we  expect  the 
theoretically  predicted  flow  speeds  for  the  developing  flow  to  be  on  the  order  of  or  somewhat 
lower  than  these  values,  which  we  seen  in  Table  V.ii  to  be  the  case.  It  is  thus  reasonable  that  the 
experimentally  observed  wheezing  sounds  were  in  fact  the  result  of  flutter  in  the  lung  airways. 


CHAPTER  7:  CONCLUSIONS  AND  DISCUSSION 


§7.1  Summary 

Motivated  by  the  fluid  dynamic  flutter  instability  seen  in  flexible  tubes  and  wheezing  lung 
sounds,  we  have  investigated  the  linear  stability  of  a  compliant  channel  conveying  a  developing 
flow.  This  flow  profile  was  approximated  by  a  plug  flow,  for  which  analytical  solutions  are 
possible,  and  with  a  perturbation  solution  (valid  in  the  upstream  region  of  the  channel),  based  on 
the  Blasius  profile,  for  which  we  obtained  analytical  long  wave  results  and  numerical  results  for 
arbitrary  wavelengths. 

The  study  of  the  perturbation  solution  showed  the  existence  of  both  a  finite  wavelength 
(flutter)  and  long  wave  (collapse)  instability,  which  may  appear  independently  or  simultaneously, 
the  latter  case  resulting  in  a  codimension  two  bifurcation  point  from  the  base  flow.  The  long 
wave  instability  is  not  present  for  our  plug  flow  solution  ,  as  it  requires  a  non-uniform  base  flow 
profile;  for  this  reason  it  has  not  been  seen  in  previous  channel  studies.  We  expect  that  the  long 
wave  instability  corresponds  to  the  tube  collapse  seen  prior  to  oscillatory  instabilities  in  tube 
experiments,  and  hence  to  flow  limitation  in  tubes  and  the  lung. 

We  have  examined  the  effect  of  varying  different  parameters  on  the  stability  of  the  system 
and  found,  in  particular,  that  effects  destabilizing  the  flutter  instability  are  those  which  would  be 
expected  to  characterize  the  lungs  of  individuals  who  are  more  prone  to  wheezing,  thus  providing 
further  corroboration  for  the  theory  that  wheezing  is  symptomatic  of  airway  flutter.  Further,  the 
long  wave  (collapse)  instability  is  more  likely  to  occur  in  wider  channels,  and  is  hence  likely  to 
appear  before  the  flutter  instability,  as  seen  in  experiments. 

We  thus  conclude,  in  light  of  the  good  agreement  between  theoretically  predicted  and 


94 


95 


experimentally  observed  flutter  frequencies  seen  in  Chapter  6,  that  we  have  further  theoretical 
support  for  the  causal  relationship  between  airway  flutter  and  wheezing  lung  sounds.  We  have 
also  found  qualitative  support  for  collapse  and  flow  limitation  being  the  manifestation  of  a  long¬ 
wave  wall-fluid  instability.  In  the  remainder  of  this  chapter  we  discuss  some  issues  related  to  this 
work.  In  §7.2,  we  discuss  the  Tollmien-Schlichting  (viscous  fluid)  instability,  and  in  §7.3  the 
issue  of  nonparallelism  in  the  base  flow  is  considered.  We  conclude  in  §7.4  with  a  brief  outline 
of  some  possible  extensions  of  the  present  work. 

§7.2  Discussion  of  the  Tollmien  Schlichting  instability 

It  is  clear  from  the  investigation  in  §6.2  of  the  TSI  in  the  context  of  the  tube  experiments 
of  Gavriely  et.al.  (1989)  that  the  TSI  is  not  likely  to  appear  for  the  thick  walled,  air-conveying 
tubes  with  which  we  are  concerned  for  the  present  study.  However,  in  other  systems,  e.g.  when 
the  ratio  of  wall  and  fluid  densities  are  similar  (so  that  the  mass  ratio  is  closer  to  unity),  or  the 
wall  elastance  is  higher,  the  TSI  may  be  significant.  While  this  is  not  the  case  in  the  lung 
airways,  it  is  for  many  engineering  applications  and  for  physiological  applications  such  as  blood 
flow.  In  the  latter,  however,  the  flows  are  likely  to  be  fully  developed,  and  in  arterial  flows 
significantly  pulsatile.  As  noted  in  Chapter  2,  the  effect  of  compliant  boundaries  on  the  TSI  for 
plane  Poiseuille  flow  has  been  considered  by  Hains  and  Price  (1962)  and  others,  and  nonlinear 
stability  calculations  have  been  done  by  Rotenberry  and  Saffman  (1990)  and  Rotenberry  (1992). 
However,  the  stability  of  pulsatile  flows  with  compliant  boundaries  has  to  date  not  been 


considered. 
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§7.3  The  effect  of  nonparallelism  in  the  base  flow 

The  derivation  of  the  Orr-Sommerfeld  equation,  as  indicated  in  Chapter  3,  requires  that 
the  base  profile  considered  be  parallel,  that  is,  a  function  only  of  the  transverse  coordinate.  When 
we  approximate  the  developing  profile  with  a  plug  flow  the  nonparallelism  in  the  base  flow  is 
moot,  but  this  is  not  the  case  when  a  true  developing  profile  such  as  our  perturbation  solution  is 
considered.  In  this  case  we  assume  that  the  flow  is  locally  parallel,  so  that  transverse  changes 
in  the  flow  are  small  and  may  be  neglected.  This  requires  that  we  choose  the  axial  position  x<) 
at  which  the  profile  is  evaluated  so  that  the  change  of  the  boundary  layer  thickness  over  a 
disturbance  wavelength  is  small  by  comparison  to  the  channel  width.  More  precisely,  it  should 
be  the  case  that  Sfxo  +  Ax)  -  5(xo)«  1,  where  fkX,)  is  the  boundary  layer  width 
(nondimensionalized  on  half  channel  width)  a  nondimensional  distance  Xo  downstream  from  the 
mouth  of  the  channel  and  Ax  =  2rt/k  is  the  nondimensional  disturbance  wavelength.  Similarity 
analysis  of  the  boundary  layer  equations  gives  6(xo)  -  5  (v*  xf  /  fl  b*2)1'2  (see,  for  example, 
Schlichting  1955);  with  this  the  condition  for  slow  variation  of  5(Xo)  reduces  in  a  straightforward 
manner  to 


2JL. «  _L(i  +  io(_^_)1/2).  (7.D 

kRwS  25  RWS 

The  requirement  for  the  developing  flow  profile  itself  to  be  valid,  i.e.  that  the  axial  position  being 
considered  be  downstream  of  the  leading  edge  of  the  channel  wall,  is  Xo  »  (Rw  S)'1  (Batchelor 
1967),  which  is  in  general  less  stringent  than  (7.1).  For  the  developing  flow  profile  that  we  use, 
the  axial  position  Xq  considered  must  also  satisfy  the  requirement  x^  <  0.02  R*  S  (c/.  §3.3).  We 
are  not  guaranteed  in  advance,  however,  that  on  satisfying  this  condition  (7.1)  will  also  be 
satisfied,  and  in  the  long  wave  limit  (k— >0)  it  will  be  impossible  to  satisfy  (7.1).  We  therefore 
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consider  here  the  effect  of  nonparallelism  on  the  stability  calculation.  For  the  TSI  non-parallel 
effects  have  been  shown  to  be  destabilizing  for  Blasius  flow  over  a  single  wall  (Gaster  1974, 
Saric  and  Nayfeh  1975,  and  others),  and  for  developing  flow  in  a  channel  (Garg  and  Gupta  1981). 
However,  the  introduction  of  non-parallelism  also  results  in  great  sensitivity  to  the  quantity  used 
to  measure  instability  (e.g.  growth  in  axial  or  transverse  velocity,  kinetic  energy,  etc.),  so  that  it 
is  difficult  to  quantify  the  degree  to  which  this  is  the  case,  and  the  magnitude  of  the  effect  may 
be  less  than  was  previously  thought  (Fasel  and  Konzelmann  1990).  Further,  as  all  of  the  previous 
studies  have  considered  only  the  TSI,  it  remains  to  be  shown  how  non-parallelism  influences  the 
channel  wall-fluid  instabilities.  However,  the  previous  studies  suggest  that  the  variation  of  the 
stability  boundary  due  to  the  degree  on  non-parallelism  present  in  the  systems  we  consider  may 
be  fairly  small,  and  a  detailed  numerical  study  for  Blasius  flow  demonstrated  that  parallel  flow 
results  provide  a  good  leading  order  approximation  for  the  desired  stability  results  (Fasel  and 
Konzelmann  1990).  We  have  therefore  chosen  not  to  undertake  an  analysis  of  the  nonparallel 
problem  here. 

§7.5  Directions  for  future  work 

There  are  several  extensions  to  the  present  work  suggested  by  the  discussions  above  and 
in  Chapter  6.  The  influence  of  nonparallelism  in  both  the  base  flow  and  collapsed  tube  geometry 
on  the  comparison  between  the  theoretical  and  experimental  results  was  seen  in  §6.2.  It  would 
thus  be  useful  to  include  the  effects  of  nonparallelism  in  the  stability  analysis  of  the  channel  flow 
that  we  considered.  By  including  these  effects  in  a  multi-scale  analysis  (as  Saric  and  Nayfeh 
1975,  for  example),  the  present  work  would  serve  as  the  leading  order  result  in  a  perturbative 
solution  in  small  parameters  measuring  the  degree  of  contraction  of  the  tube  and  the 
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nonparallelism  in  the  flow. 

A  more  complete  analysis  of  the  developing  flow  field  going  into  the  collapsed  tube 
section,  by  extending  our  FIDAP  analysis,  for  example,  would  similarly  clarify  comparisons 
between  theory  and  experiments.  Such  an  investigation  would  show  not  only  the  nature  of  the 
actual  flow  profile  in  the  collapsed  channel  section,  but  also  the  degree  to  which  three 
dimensionality  might  enter  into  the  flow. 

The  present  stability  code  could  readily  be  applied  to  investigate  the  stability  of  other 
physiological  flows,  such  as  blood  flow.  As  noted  above,  however,  arterial  flows  are  significantly 
pulsatile,  suggesting  that  the  addition  of  an  oscillatory  component  to  the  flow  field  would  be 
physically  relevant.  For  a  small  amplitude  oscillatory  component  this  could  be  incorporated  as 
a  perturbation  to  the  present  work. 

Finally,  we  have  considered  only  linear  stability  theory.  To  investigate  the  actual 
evolution  of  the  instability  a  nonlinear  analysis  is  required.  This  might  be  approached  using  a 
weakly  nonlinear  theory  similar  to  the  investigations  of  Grotberg  and  Reiss  (1984, 1985),  or  could 
alternatively  us  a  strongly  nonlinear  long  wave  analysis  (Atherton  and  Homsy  1976).  In  the 
applications  investigated  in  Chapter  6  even  the  fluner  instability  was  found  for  fairly  long 
wavelengths  (k=0.1),  suggesting  that  a  nonlinear  long  wave  analysis  may  be  useful  in  determining 
the  nature  of  the  transition  between  the  long-wave  and  flutter  instabilities. 
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APPENDIX  A:  FINITE  DIFFERENCE  SOLUTION  FOR 
DEVELOPING  CHANNEL  FLOW 


§A.l  Problem  formulation  and  uniform  grid  difference  scheme 

In  this  appendix  we  formulate  a  finite  difference  solution  for  the  flow  in  a  rigid  channel 
that  develops  from  an  initially  uniform  profile.  In  the  boundary  layer  at  the  walls  of  the  channel 
the  flow  is  governed  by  Prandtl’s  boundary  layer  equations.  We  assume  that  the  simplification 
involved  in  the  derivation  of  these  equations,  namely  the  assumption  that  axial  diffusion  is 
negligible,  may  be  applied  throughout  the  channel  and  thus  that  they  may  be  solved  to  obtain  a 
complete  description  of  the  flow.  This  reduction  is  significant,  as  it  reduces  the  elliptic  Navier- 
Stokes  equations  to  a  parabolic  form.  We  may  then  solve  the  boundary  layer  equations  with  a 
finite  difference  formulation,  treating  the  axial  coordinate  as  a  time  like  variable  and  marching 
downstream  from  the  initial  condition  at  the  channel  mouth. 

The  boundary  layer  equations  are,  scaled  on  the  initial  (uniform)  flow  speed  S*  and  half 
channel  width  b\ 

uu  +  vu  =  -^E.  +  _L  u  and  u  +  v  =  0,  (A.l) 

'  dx  Re  yy  y 

where  y  is  the  transverse  coordinate  (y=0  at  the  wall  and  y=l  at  the  midline).  Re  a  S*  b*  /  v‘,  and 
u  and  v  are  the  axial  and  transverse  velocities,  respectively.  In  an  unbounded  domain  the  value 
of  dp/dx  is  given  by  the  inviscid  solution  for  the  outer  region  of  the  flow;  in  the  case  of  a 
channel,  however,  the  central  (inviscid)  core  is  accelerated  as  a  consequence  of  conservation  of 
mass.  Thus  the  fluid  pressure,  as  well  as  the  fluid  velocities,  must  be  found.  To  accomplish  this, 
we  impose  a  constraint  requiring  that  the  flow  rate  in  the  channel  be  constant,  namely 
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J  u  dy  =  1 .  (A.2) 

o 

The  boundary  conditions  on  the  flow  are  no  slip  at  the  wall  and  a  symmetry  condition  at  the 
midline:  u=v=0  at  the  wall,  3u/3y=v=0  at  the  midline. 

To  solve  (A.l)  and  (A.2)  subject  to  the  given  boundary  conditions,  we  introduce  a  grid 
with  axial  and  transverse  indices  n  and  j,  respectively,  where  1  <  n  <  ntot  and  1  ^  j  £  jtot.  n=l 
corresponds  to  the  channel  mouth,  j=l  the  channel  wall  and  j=jtot  the  channel  midline.  We 
replace  y  derivatives  in  the  momentum  equation  ((A.l),  first  equation)  with  the  second  order 
accurate  central  difference  formulas 
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A*1 

UJ 


a2 

ay2 


fl  +  1 


1  /  n  *1 

■  237<u'“ 

-  J-«V 

Ay 


«;.V)  +  0(Ay2) 

2  u*"  +  i£V)  +  0(  Ay2). 


(A.3) 


In  (A.3)  and  equations  that  follow  we  denote  the  value  of  variables  at  the  (nj)th  grid  point  with 


a  superscripted  n  and  subscripted  j;  Ay  and  Ax  indicate  the  grid  spacing  in  the  transverse  or  axial 


directions,  respectively.  For  the  x  derivatives  in  the  momentum  equation  we  use  a  second  order 
upwind  difference  formula. 


=  -2-(— s"+1  -  2s"  +  is"-1)  +  0(Ax2),  s  =  u.  or  p.  (A.4) 
dx  Ax  2  2  1 

Finally,  to  maintain  linearity  in  the  momentum  equation  we  introduce  second  order  accurate 

approximations  for  the  nonlinear  factors  of  u  and  v  in  the  inertial  terms  of  (A.l), 


=  2  sja  - 


s  -  u  or  v , 


(A.5) 


The  continuity  equation  ((A.l),  second  equation)  is  evaluated  at  the  j+V2th  transverse  grid 
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point,  so  that  the  second  order  accurate  central  difference  approximation  for  the  y  derivative  of 
v  is 


jLv"4,1  =  _L(v"V  -  v;*1)  +  0(Ay2).  (A.6) 

dy  j*7  Ay 

The  x  derivative  of  u  in  (A. lb)  is  evaluated  using  the  upwind  difference  formula  (A.4),  with  the 
values  of  u  at  the  j+1/2th  grid  point  being  expressed  as  the  average  of  the  values  at  the  jth  and 
j+lst  points,  so  that 
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H  \ 

+  Uj  )  + 
0(Ax2,Ay2) . 


(A.7) 


The  integral  flux  condition  is  evaluated  using  a  fourth  order  accurate  open  integration 
formula, 
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(A.8) 


Using  (A.3),  (A.4),  and  (A.5)  in  the  momentum  equation,  (A.6)  and  (A.7)  in  the 
continuity  equation,  and  (A.8)  in  the  integral  flux  condition  (A.2),  we  obtain  the  difference  system 


108 


[-2-^1  -  Re^L(2v*  -  vj-l)]u*_\l  +  [4-£l.  +  3Re(2Ujn  - 
Ay2  A y  Ay2 

u;-')]«;41  +  [-2^1  +  Re^L( 2v/  -  v;*1)]^;/  +  3 Rep"*1 
1  Ay2  Ay 

=  Re(4pn  -  p"*1)  +  /?e(2a/  -  uJ‘l)(4ujm  -  m,"'1), 


(A.9) 


with  continuity. 


AX  /  (I  +  l  H  +  l  \  lr/  *l  +  l  1**1  \ 

-^(Vy.l  -  Vj  )  =  --[(My.!  +  My  ) 

\  1  /  ^"1  #1  ~  1  \  1 
2  ( My  4 !  +  My  )  +  -  (My.!  +  My  )]  . 


and  the  integral  flux  condition 


(A.  10) 


jtot  - 1 

£  My  My*  =  yror,  (A11> 

where  the  a,  are  given  in  (A.8).  To  evaluate  the  no  flux  condition  at  the  midline,  we  note  that 
from  the  central  difference  formula  (A.  3)  the  condition  9ujtot  /3y=0  is  equivalent  to  requiring 
ujtot+i=ujlot.,.  Thus  to  implement  the  no  flux  condition  we  apply  the  momentum  equation  (A.9)  for 
the  midline  grid  point  and  use  the  given  relationship  between  ujtot.i  and  ujlol_i  to  eliminate  the 
value  of  ujtot+1. 

Equations  (A.9)  and  (A.  10)  with  boundary  conditions  Ui"=0,  the  no  flux  condition,  and 
initial  conditions  uj1=l,  p=p0  form  a  problem  for  Ujn+1  and  p"+1  of  the  form  Mt  q,  =  R„  where  the 
matrix  M,  is  tridiagonal  and  doubly  bordered,  and  the  vectors  and  R,  are  given  by 
qt  =  (u^1  •••  Uj,*1*1  pn+1)T  and  R|  =  the  right  hand  sides  of  (A.9)  and  (A.  10).  This  may  be 
solved  for  u  and  p  at  the  n+lst  axial  grid  point,  after  which  (A.ll)  with  boundary  condition 
Vjn+1=0  and  initial  condition  Vj!=0  may  be  easily  solved  for  Vjn+1.  The  problem  for  Vj"+1  is  of  the 
same  form  as  that  for  u  and  p,  M2  qj  =  R2,  but  in  this  case  the  matrix  M2  is  constant  for  all  axial 
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positions  in  the  channel,  so  that  the  matrix  inversion  need  be  done  only  once  at  the  beginning  of 
the  calculation. 

For  the  calculation  of  the  velocity  and  pressure  at  the  first  point  downstream  of  the 
channel  inlet  we  do  not  have  available  the  two  levels  of  data  for  u  and  p  required  by  the  upwind 
difference  formula  (A.4),  and  so  use  a  first  order  difference  approximation  for  this  step  only, 

JLsntl  =  _L(sb+1  -  sn )  +  0( Ax),  s  =  «  or  p.  (A.12) 

dx  Ax 

§A.2  Modifications  for  a  non-uniform  grid 

The  fundamental  difficulty  in  attempting  to  solve  the  system  outlined  in  §A.l  lies  in 
resolving  the  transition  between  the  initial  condition  that  fails  to  satisfy  the  no  slip  condition  at 
the  wall  and  the  flow  in  the  remainder  of  the  channel,  for  which  we  require  the  condition  be 
satisfied.  To  obtain  an  accurate  description  of  the  flow  we  must  thus  have  a  large  number  of  grid 
points  in  the  region  over  which  this  transition  is  effected,  namely,  the  boundary  layer  near  the 
wall.  We  do  not,  however,  need  such  a  concentration  of  points  in  the  middle  of  the  channel, 
where  the  flow  is  changing  only  slightly.  This  is  accomplished  through  the  introduction  of  a  non- 
uniform  transverse  grid  to  concentrate  points  near  the  wall.  To  do  this,  we  introduce  a  new 
variable  Y=g(y),  where  g(y)  is  defined  so  that  when  the  problem  is  solved  on  a  uniform  set  of 
points  Yj  the  corresponding  points  yJ=g1(Yj)  have  the  desired  spacing.  We  use  the  transformation 
g( y)  =  Py  +  (1  -  P)(l  -  tanh(Q( l-y))/tanh(Q))  (Fletcher  1991);  P  and  Q  are  free  parameters  chosen 
to  alter  the  degree  of  non-uniformity. 

The  introduction  of  this  new  variable  changes  the  boundary  layer  equations  (A.l)  and 
integral  condition  (A.2)  through  the  replacement  of  d/dy  and  d2/df  with  (dY/dy)0/3Y)  and 
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[(d'Y/dy^O/aV)  +  (dY/dyfffidY2)],  respectively.  This  change  does  not  alter  the  difference 
formulas  used,  however,  as  it  involves  only  the  addition  of  known  coefficients;  solution  of  the 
problem  then  proceeds  as  indicated  in  §A.l. 

§A.3  Results 

To  obtain  meaningful  results  with  the  finite  difference  code  the  axial  step  size  Ax  must 
be  chosen  to  avoid  instability;  if  the  step  size  is  too  small  the  system  becomes  unstable  as  the 
flow  progresses  downstream.  At  the  same  time,  a  small  step  is  required  near  the  channel  inlet 
to  resolve  the  transition  from  a  plug  to  a  viscous  profile.  Thus  the  first  few  steps  are  taken  with 
a  small  value  of  Ax,  and  thereafter  a  larger  step  size  is  used. 

To  test  that  the  results  produced  by  the  code  are  accurate,  a  Poiseuille  profile  is 
introduced  as  an  initial  condition.  This  profile  remains  parabolic  at  arbitrary  positions 
downstream,  with  a  the  pressure  drop  that  is  linear  in  the  axial  coordinate,  as  it  should  be. 

The  developing  flow  pressure  drop  and  axial  velocity  at  various  positions  downstream  of 
the  inlet  are  shown  for  a  plug  initial  condition  in  figure  A.l.  Some  jaggedness  is  seen  for  small 
x  in  the  plot  for  pressure;  this  is  a  result  of  the  discontinuous  change  from  a  slip  to  a  no  slip 
profile  at  the  inlet.  The  difference  between  the  pressure  drop  for  the  plug  initial  condition  and 
that  for  a  Poiseuille  flow  is  a  useful  measure  by  which  to  compare  with  the  results  of  other 
investigations.  We  find  this  difference  to  be  q=0.341  (in  non-dimensional  pressure  units);  this 
is  consistent  with  other  solutions  of  the  problem.  (For  comparison,  Schlichting  (1934)  finds 
q=0.313;  Bodoia  and  Osterle  (1961)  have  q=0.338;  and  Brandt  and  Gillis  (1966)  give  q=0.331.) 


0.0  0.5  1.0  1.5 

Axial  Velocity 


Figure  A.  1:  Finite  difference  solution  for  developing  channel  flow.  A.1A: 

Pressure  as  a  function  of  axial  coordinate  x.  Solid  curve  gives 
developing  flow  result;  dashed,  Poiseuille  flow.  A. IB:  Developing  flow 
profiles  at  axial  positions  x=0.05,  0.35,  0.88,  2.8,  5.7,  10,  and  30.  Curve 
for  x=30  is  Poiseuille  flow. 


APPENDIX  B:  FIDAP  SOLUTION  FOR  FLOW  IN  A  COLLAPSED  TUBE 


§B.l  Solution  procedure 

To  use  FIDAP,  it  is  necessary  to  specify  the  domain  of  interest  and  the  mesh  desired  for 
the  finite  element  solution;  this  is  done  for  two  different  approximations  to  the  collapsed  tube 
cross-sectional  shape.  Due  to  the  two  axes  of  symmetry  in  the  problem,  a  solution  may  be 
obtained  for  one  quarter  of  the  entire  tube  cross-section  and  the  flow  in  the  remainder  of  the  tube 
inferred  from  this  reduced  problem.  The  collapsed  tube  shape  is  inherently  three-dimensional, 
so  that  it  is  not  possible  to  resort  to  axi-symmetry  or  two-dimensionality  to  further  simplify  the 
domain  considered,  but  for  the  purposes  of  this  calculation  we  are  willing  to  consider  the  steady 
flow  problem,  so  that  the  length  of  the  tube  section  in  which  a  solution  is  calculated  may  be  short 
with  periodic  boundary  conditions.  In  figure  B.l  the  domain  in  which  the  computation  is  carried 
out  is  shown  with  the  element  mesh  used.  No  slip  or  symmetry  conditions  are  imposed  on  the 
appropriate  non-periodic  boundaries,  and  the  full  steady  Navier-Stokes  equations  solved  subject 
to  these  boundary  conditions.  To  ensure  that  the  mesh  being  used  gives  an  accurate  solution  to 
the  problem,  a  mesh  that  was  half  again  as  dense  was  used  for  comparison;  little  difference 
between  the  solutions  obtained  with  the  original  and  refined  mesh  was  observed. 

§B.2  Results 

Vector  plots  showing  the  flow  profiles  generated  by  FIDAP  are  shown  in  figures  B.2A 
and  B.2B  for  the  domains  considered.  Contour  plots  of  the  same  solutions  are  shown  in  figure 
B.3A  and  B.3B. 
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Figure  B.l  Meshed  domain  for  FIDAP  computations 
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B.2A 


B.2B 


Figure  B.2  Vector  plots  of  steady  velocity  in  collapsed  tube  cross-section.  B.2A 
and  B.2B  show  velocities  for  each  of  the  two  tube  shapes  considered. 
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